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Abstract 

For scalar, electromagnetic, or gravitational wave propagation on a background 
Schwarzschild blackhole, we describe the exact nonlocal radiation outer bound- 
ary conditions (robc) appropriate for a spherical outer boundary of finite radius 
enclosing the blackhole. Derivation of the robc is based on Laplace and spherical- 
harmonic transformation of the Regge-Wheeler equation, the PDE governing the 
wave propagation, with the resulting radial ode an incarnation of the confluent 
Heun equation. For a given angular integer I the robc feature integral convolution 
between a time-domain radiation boundary kernel (tdrk) and each of the corre- 
sponding 2/ -I- 1 spherical-harmonic modes of the radiating wave field. The tdrk is 
the inverse Laplace transform of a frequency-domain radiation kernel (fdrk) which 
is essentially the logarithmic derivative of the asymptotically outgoing solution to 
the radial ode. We numerically implement the robc via a rapid algorithm in- 
volving approximation of the fdrk by a rational function. Such an approximation 
is tailored to have relative error e uniformly along the axis of imaginary Laplace 
frequency. Theoretically, e is also a long-time bound on the relative convolution 
error. Via study of one-dimensional radial evolutions, we demonstrate that the 
ROBC capture the phenomena of quasinormal ringing and decay tails. Moreover, 
carrying out a numerical experiment in which a wave packet strikes the boundary 
at an angle, we find that the robc yield accurate results in a three-dimensional 
setting. Our work is a partial generalization to Schwarzschild wave propagation and 
Heun functions of the methods developed for flatspace wave propagation and Bessel 
functions by Alpert, Grcengard, and Hagstrom (agh), save for one key difference. 
Whereas AGH had the usual armamentarium of analytical results (asymptotics, or- 
der recursion relations, bispectrality) for Bessel functions at their disposal, what we 
need to know about Heun functions must be gathered numerically as relatively less 
is known about them. Therefore, unlike AGH, we are unable to offer an asymptotic 
analysis of our rapid implementation. 
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0. Introduction 

0.1. Background. Consider the Cauchy problem^ for the scalar wave equation, 
(1) -dfU + AU^O, 

on [to, tp]xE,^, the Cartesian product of a closed time interval and Euclidean three- 
space. First, we specify suitable initial-value or canonical data U\to and dtU\to on 
at the initial time to. Next, using the rule (0, we evolve the data until the 
final time tp, along the way generating the solution U throughout the temporally 
bounded but spatially unbounded domain [to, tp]xE^. Provided physically reason- 
able initial data, this problem is well-posed; however, it is not the evolution problem 
one typically encounters in numerical wave simulation. Usually the numerical mesh 
covers only a finite portion of . 

With the finiteness of numerical meshes in mind, consider the following more 
realistic evolution problem. Let S C E^ be a round, solid, three-dimensional ball 
determined by r < rg, with a fixed outer radius, on which we specify compactly 
supported initial data U[to and dtU[to at t = to- Again, the goal is to evolve the 
data, although now generating the solution U on the finite domain = [ip, ^f] x S 
depicted in FiG. ^ Respectively, let Eq and S^? denote the ball E at i = and 
t = tp. One element of the boundary dAi = Eq U E^ U is a timelike three- 
dimensional cylinder determined by < t < tp and r — rs- Note that 
is the history in time of the spherical spatial boundary B = 9E. As it stands, 
such an evolution problem is not well-posed, since M is larger than the future 
domain of dependence of Eq. Indeed, U[t„ and dtU[to are free data, and we have no 
control over data on E"^/Eo, the region exterior to the initial E ball. Data on this 
exterior region may contain so-called ingoing radiation which will impinge upon ^B 
at later times, affecting the solution U within Ai. Most often in numerical wave 
simulation the goal is to forbid such ingoing radiation by the choice of radiation 
boundary conditions, that is explicit rules governing the behavior of U and dtU on 
^B. Often referred to as nonreflecting boundary conditions (nrbc), for the described 
problem such conditions ideally specify that the spherical boundary B is completely 
transparent. Due to the free nature of the initial data, exact nrbc are inherently 
nonlocal in both space and time. With nrbc specified along ^B, we may refer to 
such an evolution as a mixed Cauchy -boundary value problem. 

More generally, one might consider radiation boundary conditions associated 
with some other PDE and/or different type of B boundary, say cubical or irregularly 
shaped. Refs. ^ through [201 pertain either to the described spherical problem or 
to more general radiation boundary conditions. This is certainly not an exhaustive 
list, and we point the reader to review articles 11111171^0] for more comprehensive 
listings. Although we do not attempt an extensive literature review, we make 
mention of a few approaches to radiation boundary conditions in order to put our 
work in some context. Two pioneering early works are those of Engquist and Majda 
E] and Bayliss and Turkel [H] . Each develops a hierarchy of local differential 
conditions of increasing complexity. Engquist and Majda's work is based on exact 
radiation boundary conditions as expressed within the theory of pseudo-differential 
operators, and their approach is not necessarily tied to a spherical geometry nor to 

^Wc use Cauchy problem in lieu of initial-value problem in order to reserve the latter for 
the process of generating initial data, one that requires the solution of elliptic PDE for theories 
involving constraints such as general relativity or fluid flow. 



Figure 1. Finite spacetime domain c with boundary 
dM. = Eq U U '^B. Respectively, Eg and Y,f denote the solid 
round ball E at the initial time to and the final time tp- Radiation 
boundary conditions are given on the three-dimensional timelike 
cylinder ^B. Our geometric perspective on the "quasilocal" space- 
time region M comes from Rcfs. [2] and 

the ordinary wave equation. Also considering more than just our problem above, 
Bayliss and Turkel base their approach on asymptotic expansions about infinity, for 
example the standard multipole expansion for a radiating field obeying the ordinary 
wave equation Another approach to radiation boundary conditions relies on 
the introduction of absorbing layers, and as an example we mention Ref. |10| . In the 
introduction of his review article |17j , Hagstrom describes the main advances made 
in the 1990s on several fronts related to radiation boundary conditions: (i) improved 
implementations of hierarchies, such as the ones mentioned, (ii) new absorbing 
layer techniques exhibiting reflectionlcss interfaces, and (iii) efficient algorithms for 
evaluation of exact nonlocal boundary operators. Results 1151 [TB] of Grote 
and Keller fall within this first category. An advance on the second front was 
the introduction of perfectly matched layers by Berenger |18| , while a key advance 
on the third was the rapid implementation of NRBC for spherical boundaries by 
Alpert, Greengard, and Hagstrom See also related work by Sofronov [T^ . 

Hagstrom discusses the state of the art for both fronts (ii) and (iii) in his second 
review article ^5] . 

0.2. Alpert, Greengard, Hagstrom nonreflecting boundary conditions. 

Since our investigation will follow the approach of Alport, Greengard, and Hag- 
strom (agh) which belongs to category (iii) in the last paragraph, let us describe 
it in more detail. With C//,„(i,r) denoting a single spherical harmonic mode of the 
wave field (I and m are the standard orbital and azimuthal integers), we consider 
the reduced ordinary 3 + 1 wave equation. 



(2) 



Ui 2dUi l{l + l) 



Ui=0, 



dt^ 



Qj.2 



r dr r^ 
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stemming from the We have dropped the index m on Ui since it does not 
appear in the pde. Let Ui{s,r) = C[Ui]{s,r) denote the Laplace transform of the 
mode. Provided that we assume both compactly supported initial data and a large 
enough radius r, the transform Ui{s,r) obeys a homogeneous ode, 

(3) -^0 + 2.^-[-^ + ^(' + l)]t>.=O, 

known as the modified spherical Bessel equation. Here z = sr is shorthand for 
the product of Laplace frequency s and radius r. In terms of the standard half- 
integer order MacDonald function \21\ Ki^i/2{z), we define an associated function 
Wi{z) = (772/2)^/^ exp(2;)/r;+i/2(z) more closely related to a confluent hypergeo- 
metric function. The function exp(— is the outgoing solution to 
Here Wi{z) has been chosen so that Wi(z) ~ 1 as z ^ 00, that is "normalized at 
infinity." Moreover, it turns out that Wi{z) is a polynomial of degree I in inverse 
z. For example, Wo[z) = 1, Wi{z) = 1 + z'^, W2{z) = 1 + Sz^^ + iz'"^ . 

AGH introduce a nonreflecting boundary kernel, here called a time-domain radi- 
ation kernel (tdrk) as follows. Again assuming compactly supported initial data 
and a large enough radial value, the Laplace transform of the mode Ui(t,r) takes 
the form 

(4) Ui (s, r) = ai {s)r~^ exp(-sr) (sr) , 

where a/(s) is an analytic function of s depending on the details of the initial data. 
Differentiation of the last equation with respect to r leads to the identity (the prime 
denotes differentiation in argument) 

(5) sUi{s,r) + drUiis,r) +r~^Uiis,r) ^ r-^[srW[{sr)/Wi{sr)]Ui{s,r) . 

This identity is of course valid at the fixed outer boundary radius rs , provided that 
the outer two-sphere boundary B lies beyond the support of the initial data. A 
well-known classical result, Wi{z) has I simple zeros {ki n n ~ 1, - ■ ■ ,1} which lie 
in the lefthalf plane. The ki^n are also the zeros of fi'/+i/2(z), as discussed in [T^ 
and below. As a result, the object a)/(s;rB) = srBWl{srB)/Wi{srB), a frequency- 
domain radiation kernel (fdrk), is a sum of I simple poles in the lefthalf s-plane, 

(6) w,(s;rB)=> — - — . 

, s - ki,n rB 

n— 1 ■ ' 

Whence its inverse Laplace transform u!i{t;rB), the tdrk of AGH, is a corresponding 
sum of exponentials, 

I 

(7) uJi{t;rB) = cxp (fci^nt/rs) . 

n=l 

We find the following for the first three such kernels: 

(8) u;oit;rB) = 0, 

UJi{t\rB) = -r cxp{ - t/rB) , 

uJ2{t;rB) = -rg^cxp{ - ^t / r b) [3 cos {^Vst/rB) + Vs sin {^Vst / r b)] ■ 
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By elementary properties of the Laplace transform, including the Laplace convolu- 
tion theorem, inverse transformation of the identity © yields the following robc 
(a condition of complete transparency):^ 

= rg' [ LOi{t-t';rB)Uiit',rB)dt' . 
Jo 

From a numerical standpoint, this boundary condition becomes expensive for high- 
order modes, since uJi{t]rB) is made up of I exponential factors. However, elemen- 
tary exponential identities do afford an efficient recursive evaluation of the convo- 
lution jl8| . For wave propagation on flat 2-1-1 spacetime, the analogous cylindrical 
TDRK for a given angular mode is again a sum of exponentials, but now the sum is 
over a set including a continuous sector |18| . Remarkably, this is a feature shared 
by blackhole kernels. Such continuous distributions are primarily relevant for the 
low-order modes and quite expensive to evaluate. 

AGH describe an algorithm for kernel compression, by which we mean approxi- 
mation of the FDRK by a rational function, where the approximation is of specified 
relative error e uniformly for Res > 0. The resulting rational function, itself a sum 
of simple poles, typically involves far fewer terms. As a result, in the AGH approach 
the cost of updating the numerical solution at the outer boundary B is minimized 
subject to the prescribed tolerance e (which also turns out to be a relevant error 
measure in the time-domain). Therefore, we may describe their implementation of 
ROBC for the ordinary wave equation as rapid. AGH have given a formidable asymp- 
totic analysis of their rapid implementation, proving in particular that the number 
d of poles needed to approximate the fork to within the specified tolerance scales 
like 

(10) d - 0( log ly log(l/e) + log^ + log^(l/e)) 

d£,v = 1 + 1/2-^ Qo and e ^ O"*" ^Hj . Increased performance as oo is also 
seen for rational approximation of the cylindrical kernels relevant for 2-f 1 wave 
propagation. However, as remarked upon later in Section I3.3.1I for both the 24-1 
and blackhole scenarios rational approximation of low-order kernels (associated 
with costly continuous sectors) also leads to savings. 

0.3. Problem statement. Now consider the evolution problem for the scalar wave 
equation. 

describing a field U propagating on a Schwarzschild blackhole determined by g^i, 
metric functions. A slight modification of yields a wave equation flexible 
enough to also describe propagation of electromagnetic or gravitational waves on 
a Schwarzschild blackhole )23[ 1241 1^ I26| . As a model of gravitational wave prop- 
agation, the problem has applications in relativistic astrophysics: non-spherical 
gravitational collapse and stellar perturbations, among others. Gravitational wave 
propagation is also of considerable theoretical interest in general relativity. With 

'^Domain reduction appears in the early wovk |7| of Gustafsson and Kreiss, although domain 
reduction via Laplace convolution appears shortly thereafter in the work |5| of Hagstrom. In 
Ref. 1221 Friedlander considered essentially the same convolution kernel but in a different context. 



(9) 



dt 



m 

dr 



Ui 
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numerical wave simulation on a finite mesh in mind, we again choose a finite do- 
main S, now a round, three-dimensional, thick shell also bounded internally by the 
blackhole horizon H. The outer boundary B, one element of 9S = £? U iJ, is again 
specified by r = r^, while the inner boundary H corresponds to r ~ 2m (twice 
the geometrical mass of the blackhole). Let us set the task of evolving data U\tQ 
and dtU\to given on Eqi in order to generate the solution U on the finite domain 
M = [io, tp] X E with boundary dM = Eq U Sf U '"^B U ^i?.^ Here ^H, a portion of 
the future event horizon^ is the three-dimensional characteristic history [tQ,tp\y. H 
of H. To accomplish the task, we need explicit outer boundary conditions on 

= [to^tp] X B, ones stemming from the assumption of trivial initial data on the 
outer spatial region exterior to Sq- We refer to these as radiation outer boundary 
conditions (robc).^ The ROBC corresponding to Ullfl are more subtle than simple 
nonreflection, in part due to the back-scattering of waves off of curvature. 

In discussing finite outer boundary conditions, and ROBC in particular, for rel- 
ativity, we should first make a distinction between general relativity, in which the 
dynamics of spacctime is governed by the full nonlinear Einstein equations, and its 
perturbation theory, in which the dynamics of disturbances on a fixed background 
solution to the Einstein equations is governed by a linear pde similar to In 
this second paradigm one examines the propagation of weak gravitational waves on 
a fixed background spacctime (which may or may not be curved). York's survey 
article [27| on the dynamics and kinematics of general relativity is the best jump- 
ing off point for a study of the literature we now mention. Within the context of 
a mixed Cauchy-boundary value problem, Friedrich and Nagy have made theoreti- 
cal progress towards solving the full Einstein equations on a bounded domain |28j ; 
however, their results do not appear suited for numerical work. For the most part, 
approaches towards numerical outer boundary conditions in the full theory have 
relied either on matching Cauchy domains to characteristic surfaces (see [23 1^ 
and references therein) or ensuring that the outer boundary is at a large enough 
distance so that perturbation theory can be brought into play (see [21 02] and 
reference therein) . In this latter approach, the relevant perturbative wave equation 
is essentially however, the corresponding exact nonlocal ROBC are not used. 

A very different approach towards theoretical and numerical outer boundary condi- 
tions has been given in . However, it would seem of limited practical use, since 
it relies on the "many-fingered" nature of time in general relativity to completely 
freeze the flow of time at the outer boundary. In terms of harmonic coordinates 
Szilagyi; Schmidt, and Winicour have theoretically and numerically studied mixed 
Cauchy-boundary value problems for the Einstein equations linearized about flat 
Minkowski spacetime [Hi]. Attempts towards implementation of ROBC in numerical 
relativity have mostly relied on improved versions of the well-known Sommcrfeld 
condition, although Novak and Bonazzola have considered more general nonrcflect- 
ing boundary conditions with relativistic applications in mind |35| . The Sommerfeld 
condition is a local boundary condition which is exact for a spherically symmetric 

■^Although it hardly needs to be noted for this simple introduction, our time variable t is closely 
related to the advanced Eddington-Finkelstein coordinate v discussed in Section 14. 1.11 

^We include the adjective "outer" in our acronym in order to distinguish between boundary 
conditions at B and those at H . As shown later, setting appropriate boundary conditions at H 
is not nearly so difficult for our problem, heuristically since H acts as a one-way membrane of 
sorts. However, for dynamical spacetimes the issue of inner boundary conditions (at "apparent 
horizons") is a difficult problem in its own right. 
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outgoing wave in flatspace. Recently, Allen, Buckmiller, Burko, and Price have dis- 
cussed the effect of such approximate boundary conditions on long-time numerical 
simulation of waves on the Schwarzschild geometry |36) (we consider one of their 
numerical experiments below). For comments on such approaches as well as other 
remarks on numerical relativity, sec the review article by Cook and Teukolsky |37| . 
To date, there seems to have been no truly systematic analysis of algorithm error 
for treatments of ROBC in numerical relativity. 

We use frequency-domain methods and gather results which resemble those used 
and found in seminal work |38| by Leaver and also in related work |39| by Anders- 
son. Despite this resemblance, neither Leaver nor Andersson considered boundary 
integral kernels belonging to the finite timelike cylinder '^B. The starting point 
for these authors was the exact solution to the Cauchy problem as expressed via 
an integral (Green's function) representation involving Eq spatial convolution with 
initial data (actually Leaver also considered more general driving source terms be- 
yond just initial data). As Leaver noted in footnote 21 of the seeds of this 
approach are found in Morse and Fcshbach's 1953 treatment 001 of the ordinary 
flatspace wave equation, although they have origins in the 19th century. Many 
authors have since used frequency-domain and complex-analytic methods to ex- 
amine the Cauchy problem for perturbations on the Schwarzschild geometry from 
this Green's function perspective, and Andersson's article [321 is a salient recent ex- 
ample. However, we stress up-front that our problem of imposing ROBC via domain 
reduction is not the same as Leaver and Andersson's problem, and our work has 
quite a different focus on '^B temporal integral convolution. Moreover, the methods 
— those of AGH — that we describe and use in this article and its follow-up were 
only fully developed for the ordinary wave equation in the late 1990s. Section 
11.4.31 further compares and contrasts our theoretical analysis with that of Leaver 
and Andersson. 

0.4. Overview of results. In this article we describe both the exact ROBC for 
(|ll|l and an algorithm for their rapid numerical implementation. As mentioned, 
our approach to the problem follows AGH quite closely. The equation (|ll|l is lin- 
ear, but necessarily with variable coefficients. Nevertheless, exploiting its time and 
rotational symmetries, we may likewise use Laplace and spherical-harmonic trans- 
formation in order to obtain a second-order radial ode which turns out to be an 
incarnation of the confluent Heun equation [411 142j , also related to the generalized 
spheroidal wave equation^ discussed in some literature |43l I44| . Following AGH, we 
may formally introduce the tdrk as the inverse Laplace transform of the homoge- 
neous logarithmic derivative of the asymptotically outgoing solution. Analytically, 
the FORK, the logarithmic derivative in question, is a sum of poles, although now 
the sum is over both a discrete set and a continuous set (similar to the situation 
for the flatspace wave equation in 2 -I- 1 rather than 3+1 dimensions \1H\). 

Employing both direct numerical construction of blackholc fork's as well as 
their compression along the lines of AGH, we numerically implement the exact ROBC 
for waves on Schwarzschild. In principle, we may implement the conditions to 
arbitrary numerical accuracy even for long-time simulations. Via study of one- 
dimensional radial evolutions, we demonstrate that the described ROBC capture 
the long-studied phenomena of both quasinormal ringing and late-time decay tails. 



°The ordinary spheroidal wave equation stems from variable separation of the ordinary wave 
equation Q in oblate or prolate spheroidal coordinates I4UI . 
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Since our implementation is based on the exact nonlocal history-dependent ROBC, 
we sidestep issues raised by Allen, Buckmiller, Burko, and Price and are indeed 
able to capture decay tails with our boundary conditions. Wc demonstrate this 
below with an example considered in |36j . We also consider a three-dimensional 
evolution based on a spectral code, one showing that the ROBC yield accurate results 
for the scenario of a wave packet striking the boundary at an angle. Our work is 
a partial generalization to Schwarzschild wave propagation and Heun functions 
of the methods developed for flatspace wave propagation and Bessel functions by 
AGH, save for one key difference. Whereas AGH had the usual armamentarium of 
analytical results (asymptotics, order recursion relations, bispcctrality) for Bessel 
functions at their disposal, what we need to know about Heun functions must be 
gathered numerically as relatively less is known about them. 

Due to this key difference, we are, unfortunately, unable to offer a rigorous 
asymptotic analysis of our rapid implementation. However, our numerical work 
suggests that the number d of approximating poles grows at a rate not at odds 
with the one mentioned above. Indeed, as seen in Table El from Section 13.3. II for 
e — 10"^" and tb = 30m we have found that d = 20 is sufficient for all I < 64. 
With the same rs but e = 10~^ instead, d = 14 is sufficient for all I < 256. These 
d values gives us some idea of the cost associated with our rapid implementation. 
Let us focus attention on a single spherical harmonic mode, ignoring the cost of 
performing a (numerical) harmonic transform. Depending on the algorithm used for 
wave simulation, such a transform may or may not be performed at each numerical 
time step. For the range of I values and e accuracies we consider, d is roughly on the 
order of 10. Therefore, at each numerical time step we expect a boundary operation 
count which is some multiple lOp of this d value. (Actually poles whose locations are 
properly complex are more costly than those which lie on the real axis. This affects 
p, but no matter.) A practical radial discretization of [2m, 30m] will have some 
multiple lOOOg of a thousand mesh points, whence p/q is a straightforward estimate 
for the percentage cost (at each numerical time-step) of updating the solution at the 
edge of the computational domain relative to the cost of updating on the interior. 
Note that while p remains fixed, q increases with mesh resolution, so that the cost 
of our ROBC relative to interior cost becomes accordingly negligible. Let us support 
this assertion with a concrete example. For the one-dimensional radial evolution 
described in Section 15.1.11 we have kept track of the total CPU time spent both 
on interior work and the ROBC. We list the ratio of these times as percentages in 
Tabled Such percentages belie the true savings of the implementation. Without 
some form of ROBC, one would be forced to consider the free evolution of waves on 
a domain larger than [2m, 30m] , one large enough to ensure that the waves would 
not reach the larger outer boundary during the simulation.^ The cost of using a 
larger domain as "boundary conditions" is usually some multiple of the interior 
cost. Therefore, our ROBC typically cost less than a percent of what evolution on 
such a larger domain costs. In SECTiON l4.3.3l we discuss memory and storage issues 
relevant to our implementation of ROBC. 

0.5. Summary. Section Q] discusses variable transformations, various resulting 
forms of Hll|) . and the exact ROBC. We start off by defining dimensionless coordi- 
nates for time r, radius p, and Laplace frequency a. For example, r = 2mp and 



^Or at least large enough so that waves reflected off the outer boundary would not disturb the 
smaller computational domain of interest. 
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I 


Number of radial mesh points 


1024 


2048 


4096 


8182 


16384 





2.24 


1.49 


1.00 


0.39 


0.21 


33 


2.41 


1.76 


1.31 


0.50 


0.26 


64 


2.55 


1.81 


1.37 


0.51 


0.27 



Table 1. Relative cost of e = 10^^" robc. For the one- 
dimensional evolution described in Section I5.1.1I we list percent- 
age values for CPU time spent on robc relative to CPU time spent 
on interior work. Whence each value is essentially the aforemen- 
tioned ratio p/q for a particular q. 



s ~ a /{2m). The outer boundary B is determined p ~ pB- With these coordinates 
we introduce the asymptotically outgoing solution Wi{ap; a) to the radial equation, 
the one corresponding to the Bessel-type function Wi{sr) = Wi{ap) above. For a 
given angular index I the tdrk uji{t;pb) is the inverse Laplace transform of the 
FORK uji{a;pB) = (TpBWI{apB',o-)/Wi{apB',cr)- We then write the ROBC as an in- 
tegral convolution between uji(t;pb) and each of the corresponding 2/ -I- 1 modes 
ipimiT, pb) of the radiating field ■0(r, p, 9, (j)), where ip isU from above but expressed 
in terms of different coordinates. Afterwards, we describe the key representation 
of 6ji{a;pB) as a (continuous and discrete) sum of poles. Section ^ ends with 
the derivation of an estimate for the relative error e associated with approximating 
'^liT', Pb) by a numerical kernel £,i(t] pb)- 

Section El describes numerical evaluation of both Wi{(TpB]cr) and u}i{a;pB), 
with the former considered as a function of complex Laplace frequency a (mostly 
lying in the lefthalf plane) and the latter as a function of purely imaginary a = iy. 
Both types of evaluation rely on numerical integration over certain paths in the 
complex plane. We consider several numerical methods, but the main ones involve 
path integration in terms of a complex variable z = a p. While numerical evaluation 
of Wi {(jpB] c) is important insofar as studying the analytic structure of this function 
is concerned, implementation of ROBC mainly requires that we are able to accurately 
evaluate 6ji{iy; pb) for any y G M.. In this section we also discuss in detail the 
accuracy of our numerical methods. 

Section 01 focuses on the sum-of-poles representations of u}i{(j; pb)- The first 
subsection is a qualitative description of the analytic structure of Wi {crpB ; f) and its 
relevance for the exact sum-of-poles representation. This subsection examines the 
zeros in frequency cr of Wi{apB;(j) which correspond to poles of uji{a;pB)- It also 
studies the branch behavior of Wi{apB', cr) along the negative Recr axis, behavior 
that gives rise to a continuous pole distribution (these are not really poles in the 
sense of complex analysis). This distribution appears in the exact sum~of-poles 
representation, and we graphically examine it. The second subsection presents our 
direct numerical construction of (jJi{a; pb) for < Z < 10 and Pb € [15,25]. We 
discuss several numerical accuracy checks of our direct construction. The third 
describes kernel compression, and the resulting approximation of ilii^a; pb) by a 
rational function which is itself a sum of poles. In this third subsection we consider 
the bandwidth < / < 64. 
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Section 0] presents the details of our implementation of robc. As mentioned, 
our implementation is designed around the MacCormack predictor-corrector algo- 
rithm )45| . The first subsection describes both our choice of spacetimc foliation 
and the associated first-order evolution system of pde, while the second reviews 
the MacCormack algorithm in the context of this system. As an aside, the second 
subsection also addresses the issue of inner boundary conditions at the horizon H . 
The third subsection describes the implementation, showing how robc fit within 
the framework of the interior prediction and correction. In this final subsection we 
also address some memory issues relevant to our implementation. 

Section |5l documents the results of several numerical tests of our implemen- 
tation of robc. Throughout this section, we consider a blackhole of mass m = 2 
enclosed within an outer boundary of radius = 60 so that the horizon is located 
at 2m = 4 and pB = 15- The choice m = 2 is not particularly special, and has 
been made only to provide an example in which the mass is neither ^ nor 1. The 
bulk of this section focuses on one-dimensional radial evolutions of single spherical- 
harmonic modes. However, in the last subsection we consider a three-dimensional 
test. 

A final sections offers some conclusions and also discusses applications and ex- 
tensions of our results. 

In Appendix we discuss a certain modification of the MacCormack algorithm 
which is featured in our implementation of robc. In Appendix ^ we provide 
numerical tables for several compressed kernels. 

Starting on page |^ we have listed our main symbols in order of appearance, 
with the section number given for where each symbol first appears. Symbols not 
listed there are defined and used locally. Although some symbols in this work have 
multiple meanings, within a given section a symbol's meaning does not change. We 
point out that as a complex variable a = x -\- iy, and, therefore, for the complex 
variable z = ap always write z = Rez -|- ilmz. 

1. Wave equation and radiation outer boundary conditions 

This section sets up the theoretical framework on which the subsequent sections 
rest. We here discuss various pde and ode relevant for wave propagation on the 
Schwarzschild geometry. We then derive the exact and nonlocal radiation outer 
boundary conditions (robc) appropriate for asymptotically outgoing fields, thereby 
paving the way for their numerical implementation in later sections. 

1.1. Wave equation on Schwarzschild background. 

1.1.1. Line-element. Consider the diagonal line-element describing a static, spher- 
ically symmetric, vacuum blackhole of mass m, 

(12) ds^ = -FdT^ + F^Mr^ + r^dd"^ + sin^ ed(f>^ , 

written with respect to the standard static time T and areal radius r [521 . We 
use uppercase T here to save lowercase t for a different time coordinate needed 
later. Note that the metric coefficient F{r) = 1 — 2m/r vanishes — and so is 
singular — as r ^ 2m. As is well known, r — 2m does not represent a physical 
singularity, rather the coordinate system is degenerate for this value of the radius. 
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In these coordinates the round sphere determined by r = 2ni is the bifurcate cross- 
section of the event horizon of the blackhole. In this work, we are chiefly interested 
in the "exterior region" defined by 2m < r < oo. 

It will prove convenient to pass to and work with dimensionless coordinates 
(r, p, 6, (j>) defined by 



(13) 



2mr = T , 2m/9 



(6 and (j> are already dimensionless). After the rescaling ds^ i—s- ds^/(4m^), we may 
rewrite the line-clement (I12II in the dimensionless form 



(14) 



ds^ = -Fdr^ + F-^dp^ + p^d9^ + p2 sin^ 



where now F{p) = 1 — 1/p so the unphysical singularity is located at p = 1. 

We also consider the outgoing and ingoing systems of Eddington-Finkclstein 
coordinates |52[ I53| . here in dimensionless form. To construct them, first introduce 
the Regge- Wheeler tortoise coordinate |24[ I52| 



(15) 



p*= p + log(p - 1) . 



Recall that this transformation is valid for 1 < p < oo which corresponds to — cx) < 
< oo, and that p — > 1+ corresponds to p* — oo. In terms of the tortoise 
coordinate we write H14() as 



(16) 



ds' 



F{-dT^+dpl 



9^d9^ 



p^ sin^ 



The characteristic coordinate p = t — p^ is the retarded time, and the set (p, p, 9, (f) 
is the outgoing Eddington-Finkelstein system. With respect to it, the line-element 
takes the form 



(17) 



ds^ = -Fdp^ - 2dpdp 



o^d9'^ 



p^ sin^ 



In the {p,p,9,(f)) system gpp — 0, so that the vector field d/dp is characteristic or 
null, whereas in the (r, p, 9, <f)) system gpp ~ F~^[p), so that d/dp is spacelike on the 
exterior region. Level-/z hypersurfaces are characteristic and outgoing (cones which 
open up towards the future) with d/dp as their outgoing generator. In Section 
14.1.11 we consider the ingoing Eddington-Finkelstein coordinate system which is 
based on the characteristic coordinate = r + p* known as advanced time. 



1.1.2. Wave equation. The covariant d'Alcmbcrtian or wave equation associated 
with the diagonal line-element H14() is the following: 



(18) 



1 - 



dr^ 



p2 dp 



1 - 



1\ (hp_ 
dp 



= 0, 



where A52 is the Laplace operator (with negative eigenvalues) belonging to the 
unit-radius round sphere 5*^. Notice that we use tjj for the wave field associated 
with the static time coordinate r (or associated with its counterpart T) introduced 
above, whereas in the introduction we have used U for the wave field. Later we 
use U for the wave field associated with a certain time variable t related to ingoing 
Eddington-Finkelstein coordinates. Our numerical work is based on t (which is 
why U and t, rather than ^p and T, appear in the introduction). For flat spacetime 
T and t are the same, and so ijj and U are also formally the same for m — 0. 



11 



Introducing the standard set Yim{0^(j)) of basis functions for square-integrable 
functions on 5^, we consider an appropriate expansion 



(19) 



of the field ip in terms of spherical-harmonic modes ij^im ■ The spherical-harmonic 
transform of (|18|l . 



(20) 



1 - - 



dp 



1 

1 - - 
P 



dp 



p2 



0, 



is the PDE governing the evolution of a generic mode "0; . On Tpi we have suppressed 
the m, since it does not appear in the pde. 

Addition of a single simple term to (|20|) yields a modified wave equation flexible 
enough to describe either the mode evolution of an electromagnetic field Ap or 
the mode evolution of small gravitational perturbations Sgap on the Schwarzschild 
background. The modified equation is 



(21) 



1-i 

p 



9r2 



P 



]_d_ 

2 dp 



P 



p) dp 



= 



with the spin j = 0, 1, 2 corresponding to scalar, electromagnetic, and gravitational 
radiation. Wc review the history of this correspondence in the next paragraph. We 
may cast ()21|l in a particularly simple form via simultaneous transformation of the 
independent and dependent variables. Indeed, setting 'I'; = prpi and here viewing 

)d/dp, we rewrite (|21|l as follows: 



d/dpf, as shorthand for (1 — p ^ 



(22) 



9r2 



dpi 



The Regge-Wheeler potential 
(23) V{p) = 



1 

1 - - 
P 



1(1 + 1) 



p- 



would depend only implicitly on were we using as the independent variable. 
As we will see in Section 11.2.31 the Laplace transform of (|22|l is important theo- 
retically, since it elucidates the role of Laplace frequency as a spectral parameter. 

Wheeler derived the j = 1 version of I|22I23|I in 1955 [221, showing that each of the 
two polarization states for an electromagnetic field on the Schwarzschild geometry 
is described by one copy of the equation. Regge and Wheeler then derived the j — 2 
equation for odd-parity (or axially) gravitational perturbations in 1957 |24| . and 
Zerilli introduced a similar equation describing even-parity (or polar) gravitational 
perturbations in 1970 |2SI- In the 1970s Chandrasekhar and Dctwciler demonstrated 
that the Zerilli equation can be derived from (|22|1 . although the derivation involves 
differential operations (see and references therein). 

Adopting p as the time coordinate, we define ipi {p, p) = ipi{p + p^,, p) and write 
(EU as 



(24) 



dpdp 



2d^ 
P dp 



P 



)_d_ 

'2 dp 



1 

1 - - 
P 



d(pi 
dp 



P' 



p3 



0. 



Another way to obtain this wave equation is to form the d'Alembertian associated 
with H17|) and then implement a spherical-harmonic transformation. Similar to 
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Static time coordinate system (r, p) 


Retarded time coordinate system (/i, p) 


i'li-r, p) 


ODE for L.t. is analogous to 
the spherical Bessel equation. 




J = ODE for L.t. is directly 
the confluent Heun equation. 


p) 


ODE for L.t. elucidates role 
of (T as a spectral parameter. 




ODE for L.t. has outgoing 
solution normalized at oo. 



Table 2. Wave fields and their relevance. L.t. stands for 
Laplace transform. As wc discuss in Section I1.2I via the Laplace 
transform we trade a pde for an ode. 



above, we may either set $; = pcpi or p) — '^i{p, + p^, p), thereby expressing 
(PI as 

(25) 2f-^-^ + np)*; = 0, 

dpdp* dpi 

again viewing djdp^ as a shorthand. As we will see, for a given / the fdrk is 
built from the outgoing solution to the formal Laplace transform of H25|) . Table El 
lists the wave fields we have introduced, and it also briefly describes the theoretical 
importance of each field's Laplace transform. The statements made in the table 
are explained in Sections 11.21 and ll. 31 

1.2. Laplace transform and radial wave equation. 

1.2.1. Laplace transform. The Schwarzschild geometry is static,^ and with respect 
to the chosen coordinates we indeed sec that the components of the metric tensor 
are r-independcnt. In turn, the variable coefficients of the linear wave equation 
described in the last subsection do not depend on time, a scenario permitting 
study of the equation via the technique of Laplace transform. The description 
of the technique in Chapter 12 of the textbook [S^l by Greenberg is suitable for our 
purposes. Let C denote the transform operation, 

/•oc 

(26) C[9]{a)= / e"^^g{T)dT. 



Jo 

Here we use cr for the variable dual to r with respect to the Laplace transform. 
We may define a physical variable s = cr/(2m), with dimensions of inverse length, 
which is dual to t and satisfies st ~ err. We may also define a formal Laplace 
transformation on the retarded time p by replacing r with p in the last equation. 
For the time being we proceed with the transformation on r. 

1.2.2. Laplace transform of the wave equation. Let us formally compute the Laplace 
transform of (|21|l . in order to get an ODE in the radial variable p. With ipi = C[ipi] 
and a dot denoting t differentiation, we have 

'C[^/'i](cr,p) = (Tipi{cT,p) - lpi{0,p) , 

(27) C[i'i]{cT,p)^a^iia,p)-aMO,p)~MO,P), 

provided Recr > 0. We assume that the initial data ipi{0, p) and '0i(O, p) vanish in a 
neighborhood of p, as is true for compactly supported data so long as we choose p 



'''More precisely, d/dr is a hypersurface-orthogonal vector field which satisfies the Killing 
equation i^a/SrS/Ji' = Oj where £ denotes Lie differentiation. 
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large enough. This assumption ensures formaUy that upon Laplaee transformation 
we may replace t partial differentiation by cr multipUcation. Whence, after some 
simple algebra we find 



{p-ir P{p-l) pHp-1) 



i'i = o 



dp2 p[p - 1) dp 

for the Laplace transform of H21|l . 

It is instructive to see what happens to 128|) in the m 0+ limit. Before taking 
the limit, first recall that p = r/(2m) and a = 2ms, so that the product z = up = sr 
is independent of m. With this in mind, we divide the overall equation by a factor 
of cr^ and find 



(29) 



d V/ 2z — cr Aipi 
dz-2 z{z — a) dz 



2 



+ 



^/ = 0, 



(z — cr)2 z(z — cr) z'^{z — a) 

where a is now shorthand for 2ms. Formally then, the m ^ 0"'' limit along with 
multiplication by sends the last equation into the modified spherical Bessel equa- 
tion (msbe) pnj : 

(30) ^'^^ + - W + l{l + 1)] V3/ = . 

dz^ dz '■ ■' 

As linearly independent solutions of the msbe we may take 

(31) ki{z) = ' *'(^) = \/^^'+i/2(^) ' 

where Ki^i/2iz), MacDonald's function, and /;+i/2(z) are standard modified (cylin- 
drical) Bessel functions of half-integer order (21]. Later we emphasize the close 
parallels between these Bessel functions and the corresponding solutions to (|29|l . 

1.2.3. Laplace frequency as a spectral parameter. We observe that 

(32) ^ - V{p)h = a'^i , 

is the formal Laplace transform of H22|) . This is a remarkable form of the radial 
ODE for several reasons. First, as might be expected from the suggestive form 
of the equation, we could consider it in the context of an eigenvalue problem, 
although one in which the operator on the lhs is not self-adjoint. More precisely, 
suppose we seek solutions to H32(l which vanish a.t p = ps (a fixed constant) and 
are also asymptotically outgoing, that is behave as e~'^'' for large p. We do then 
(numerically) find solutions corresponding to a discrete (but finite) set of a values, 
but these turn out to be values in the Icfthalf plane. ^ For such a the term e~°''' 
blows up as p gets large, spoiling any possible self-adjointness for d^/dpj — V{p) 
on [pb, oo) with these boundary conditions. 

Let us also consider the Bessel analog of (I32II . Namely, 



(33) ^ - '-^^^i = a'-^ 



^([+l).f, _2. 

8n 



The results we describe in subsequent sections justify this statement, although in what 
follows we work with a different form of the ODE stemming from yet another transformation 
^'j = exp(— (Tp«)$( of the dependent variable. See Section II. 3. 21 and what follows. 
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We reach this equation by first passing to z = ap as above, taking the m ^ 0+ hmit, 
and then passing back to p = z/a. For the type of eigenvalue problem mentioned 
above, the operator (P/dp^ — l{l + l)/p'^ is again not self-adjoint; however, this fact 
is not our prime concern now. The discussion in SECTION II . 2 . 21 shows that H33() has 
solutions, such as {c^p)^^'^ Ki^i^2i'^p)i of a special form. Indeed, they simultaneously 
solve an ode in the spectral parameter a |56j . 

(34) ^-^i. 



Accordingly, we describe solutions to H33|l as bispectral. Unfortunately, solutions 
to the more complicated ODE H32|) are not bispectral in this sense, and the lack 
of an associated differential equation in the spectral parameter a complicates our 
numerical investigations. More comments on this point follow in Section I2.1.2I 

1.3. Normal and normalized form of the radial wave equation. 

1.3.1. Normal form. Standard analysis jS3 of the ode H28() shows that p — 
and p = 1 are regular singular points, corresponding respectively to indicial 
exponents ±j and ±cr, whereas p = cx) is an irregular singular point |54| . To put 
(|28|) in a "normal form," we transform the dependent variable ■(/;/ in order to (i) set 
one indicial exponent to zero at each singular point and (ii) "peel-off" the essential 
singularity at infinity as best we can. To this end, let us set 

(35) i'l = p^ip-iye-'^PQi = p^e-^P'Qi , 

noting that such a transformation clearly achieves condition (i). Moreover, the 
large-p behavior of suggests that in order to achieve condition (ii) we should 
peel off either the factor exp(— crp) or exp(tTp). Our choice of peeling off exp(— crp) 
indicates our intention to examine asymptotically outgoing radiation fields. In H35(l 
we have peeled off a factor of {p— rather than (p— 1)"^ in order that the tortoise 
coordinate appears in the argument of the exponential factor in the transformation. 
Under our transformation Eq. (|28|l becomes 
(36) 



d^9; 

dp2 



-2cr 



1 + 2j 1 - 2cr 



P 



dp 



8; = 0. 



p-1 p(p-l) 

We remark that one may also obtain the j = version of (|36|l directly from H24() via 
formal Laplace transform on the retarded time fi, i. e. for j = we can say &i = (pi. 
Eq. H36(l is a realization of the (singly) confluent Heun equation |41[I42| 



(37) 



d^G 
dp2 



+ ^ 
P 



1 



dG 
dp 



a/3 



p-1 p(p-l) 



G = 0. 



which has the generalized Riemann scheme |42| 

2 



(38) 



1 



1-7 



1 
1 


l-<5 



;p 



7 ■ 



oo 
a 

- 5 — a 


-/3 

The first three columns of the scheme's second row indicate singular-point locations, 
while the corresponding columns of the first row indicate their types. That is to say, 
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we have regular singular points at p = and p = 1 and an irregular singular point 
at oo which arises as the confluence of two regular singular points (the 2 in the third 
column of the first row indicates this confluence).^ The remaining information in 
the first two columns specifies the indicial exponents at the regular singular points, 
while the remaining information in the third column specifics the Thome exponents 
corresponding to the two normal solutions about the point at oo. These solutions 
have the asymptotic behavior 

(39) G+ - p-'^e^-P , G- - p-<-s+°'e-'^P , 

as p oo (in some sector which we discuss later). Finally, in the fourth column 
we have the independent variable p as well as the accessory parameter q. An ode 
with the singularity structure of the confluent Heun equation is determined by 
the indicial exponents belonging to the regular singularities along with the Thome 
exponents only up to a free parameter q. Appendix B of ^ discusses this point in 
more detail. 



1.3.2. Normalized form at infinity. Viewed as the confluent Heun equation, we see 
that our radial wave equation H36|l has the following generalized Riemann scheme: 

1 1 2 

1 oo ] p 

1+j ;jij+l)-lil + l) 
-2j 2a \+]-2a 


2a 



(40) 



showing that the normal solutions to l|36|) obey 
(41) 



e+ ^ p- 



er ~ p- 



-J+2o"g2crp 



as p — !• oo. Wc may also write ^ P~ cicp{2ap^) for the large~p behavior of the 
second solution. The scheme (|38|l shows that confluent Heun functions are generally 
specified by five parameters. However, our specific scheme (|4n(l corresponds to a 
two-parameter family of functions (those parameters being a and I, with j viewed 
as fixed). This is comparable to the situation regarding the flatspace radial wave 
equation and the associated one-parameter family of Bessel functions (that param- 
eter being the Bessel order I + 1/2). Bessel functions (suitably transformed) are a 
one-parameter family within the larger two-parameter class of confluent hyperge- 
ometric functions (which may also be represented as either Whittakcr functions or 
Coulomb wave functions) (581 1591 1^ . 

Numerical considerations below dictate that wc work instead with an outgoing 
solution which is "normalized at infinity," that is to say approaches unity for large 
p. Therefore, we now enact the transformation 



(42) 

or in terms of the original field t/ji ~ p^^ exp(— (Tp*)$; 



(43) 



d^ 

dp2 



-2a- 



1 - 2a 
P-I 



d$, 
dp 



1-3' 



whereupon we find 

l-3^ + l{l + l 



=0 



^Appendix B of Q] shows how the confluent Heun equation arises from the Heun equation, an 
ODE similar to the familiar hypergeomctric equation, although possessing four rather than three 
regular singular points. 
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as the ODE satisfied by $/. Remarkably, tliis equation agrees with that obtained 
directly from (|25(l via formal Laplace transform on the retarded time whence 
our choice $/ with a hat for the dependent variable here. We emphasize that this 
statement is true for all possible spin values (j = 0, 1, 2), whereas the identification 
Qi = ipi mentioned before is valid only for j = 0. 

We again set z = ap = sr (independent of m) and divide (|4^{|l by an overall 
factor of a^, thereby reaching the following particularly useful form of the radial 
wave equation: 

1-f + 1(1 + 1)' 

dz 



d2$ 

(44) -T-I 



„ 1 1 - 2a 
-2-- + 



$i =0. 



z{z — a) 

With Wi{z] a) and Zi[z; a) respectively denoting the outgoing and ingoing solutions 
to this ODE, the corresponding solutions to H43|l are Wi{ap;a) and Zi{ap;a) with 
(T here viewed as fixed. As p — > oo these obey 

(45) Wi [ap; a) ^ 1 , Zi {ap; a) ^ e^-"* , 

as shown by the material presented above. Respectively, we might also denote Wi 
and Zi by and We will often refer to Wi{crp; a) as a confluent Heun function 
(or just Heun function), even though it differs from a Heun function by the p^^^-' 
factor (as shown above). This terminology streamlines our presentation, allowing 
us to draw the flatspace/Schwarzschild distinction via the Bessel/Heun modifiers. 

Recall that ct — *■ as m ^ 0^, whereas the product z ~ ap remains fixed in the 
said limit. Therefore, in the m ^ 0+ limit Eq. 1)44(1 becomes an ode 

(40) f*l-2«i-i(i±ll4,^o 

dz^ dz z^ 

which could also be obtained straight from the msbe ((3()|1 via the transformation 
ifji = z~^e~^$i. In terms of the two-parameter functions introduced above, Wi{z] 0) 
and Zi{z; 0) are respectively the formal outgoing and ingoing solutions to ((46(1 . We 
shall also write these as simply Wi{z) and Zi{z) when there is no cause for confusion. 
A few examples may be illuminating. The first three outgoing solutions to the msbe 
are the following spherical MacDonald functions: 

l\ , , . 3 3 



(47) fco(z) = , ki{z)^ l + k2{z)^ 1, , 2 

z z \ z / z \ z z^ 

Now consider the following polynomials in inverse z: 

(48) Wo(z)-l, Wi{z)^l + -, W2{z) = l + - + \. 

Z Z Z'^ 

From the discussion above we see that these are outgoing solutions to 1(46(1 . and 
clearly ones which are normalized at infinity. We shall see that outgoing solutions 
Wi{z\ a) to 1(44(1 are similar, albeit not simple polynomials in inverse z = a p. 

1.3.3. Asymptotic expansion for normalized form. For our purposes Eq. 1(44(1 will 
prove the most useful form of the frequency-space radial wave equation, so let us 
describe its outgoing solution Wi (z; a) as a formal asymptotic series. Our discussion 
in Section II .3.21 focused on the variable p, but the same normalization issues are 
pertinent for z. First, for convenience we set n = 1 — j^; hence k takes the values 
1,0,-3 for scalar, electromagnetic, and gravitational cases respectively. We will 
often work with k instead of j. Here and in what follows, we suppress the solution's 
K dependence. 
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Assume a solution to ()44|l taking the form 

oo 

(49) W^,(z;a)~^d„(a)z-", 

ri=0 

demanding that do{a) ~ 1. Of course the remaining d„(<T) wiU in general also 
depend on / and k, but we suppress this dependence here. Standard calculations 
then determine both di(cr) = l{l + l)/2 and the following three-term recursion 
relation: 

,^„, , , , [l{l + 1) - n{n + l)]dn{cr) + a{n^ + K - l)dn-i{<j) 

(5Uj d„+i(crj = — — — . 

2{n + 1) 

A dominant balance argument shows the dn+iic) / dn{cr) = 0{n), whence the series 
(|49|l is generally divergent and only summable in the sense of an asymptotic expan- 
sion. Olver shows that the sector of validity for this asymptotic expansion includes 
the entire z-plane (sec Chapter 7 of 161'). 

Set c„ = dn{0)- Sending a i-^ in (|r)0|l then yields the simple two-term recursion 
relation 

[l{l + l)-n{n + l)] 

51 c„+i = — — c„, 

2{n + 1) 

with solution (see Ref. [HI, P- 202) 

' ~ 2"n!r(/-n+l) ■ 

When / is an integer, as is the case here, the series X]^o truncates, showing 

for all / that the solution Wi{z) is a polynomial of degree I in inverse z. All 
coefficients c„ are positive and nonzero, and we can ultimately conclude that all 
zeros of Wi{z) lie in the lefthalf z-plane. Furthermore, the last nonzero coefficient 
is 

r(2/ + i) 

(53) c, = ^^, 

and from this formula we may appeal to the asymptotic behavior of the gamma 
function (see Ref. 001 1 P- 486) in order to show 

(54) Q~r(/)2'V^ 
as I becomes large. 

1.4. Radiation outer boundary conditions. This subsection derives and dis- 
cusses ROBC for a single spherical-harmonic mode tpim (t, p) , but we continue the 
practice of everywhere suppressing the subscript m. This subsection's formulae 
are valid for all possible spin values j = 0, 1, 2; however, for concrete examples we 
choose J ~ 0. 

1.4.1. Derivation of the radiation kernel. Although we will now derive exact equa- 
tions, let us define the radial computational domain to be 2m times the p interval 
[1,P_b]. The radial numerical mesh will be a discretization of this domain. Now 
consider an infinite radial domain S defined by p > pmax, with pmax < PB- Let So 
denote the intersection of S x [0, oo) with an initial t = Cauchy surface. Assume 
that the initial data p) ^-^d V'i(0, p) is of compact support and, moreover, van- 
ishes on Sq. The condition pmax < Pb ensures that the computational domain edge 
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PB does not intersect the support of the mitial data (see FiG.EJ. Then the analysis 
of the last subsection establishes the formal expression^'' 



(55) Mr,p)^C 



-1 



^ e-P'Wi{ap-cj) , ^ , e'-P'Zi{ap-CT) 
ai[<T) h bi(a) 



P P 

as the general solution to the wave equation (|21|l on the history S x [0,oo) of So- 
Here the coefficients ai{a) and hi{a) are arbitrary functions analytic in the righthalf 
(T-plane. Now, Zi{ap;a) cxp(2crp,) as p — > oo, showing that bi{a) must be zero 
(for otherwise the solution is not asymptotically outgoing as expected for an initial 
"wave packet" of compact support). 

The solution within the computational domain is also obtained as above via 
inverse Laplace transform. However, now the relevant frequency-space radial func- 
tion solves the inhomogcneous version of H28|l . in which case the source 

(56) {p^^pr'Jiip;a)^-p^ip-ir^[MO,p) + <yMO,p)] 

replaces zero on the rhs of the equation, as is necessary for non-trivial initial data. 
This solution appropriately matches ai{(T)p~^e~'^P'Wi{ap; a) at the largest radius 
Pmax on which the data is supported. Let us further assume that the initial data 
is supported only on [pmin , Pmax] , where pmin > 1 (again, see FiG.EJ. In this case, 
taking the second solution p"^ exp(— (Tp*)Zi(crp; cr) as ingoing at the horizon and 
using well-known methods associated with one-dimensional Green's functions, one 
can show that 

(57) ai{a) (X [ cxp{~ap){p-iy^p^^Zi{ap;a)Ji{p;a)dp, 



where the proportionality constant is determined by calculating the Wronskian of 
the two chosen linearly independent solutions to the homogeneous equation j54| . 

Let us now derive the explicit form of the radiation kernel, assuming that we now 
work in the region p > ps- Consistent with our presentation thus far, we denote 
by ipi (r, p) the function satisfying 

(58) C[^i]ia, p) ^ M<^, p) ^ ai{a)p~h"'P'Wi{ap;a) . 

Differentiation of this formula by p then gives 

Wl{crp;<T) a 1 



(59) dpM'^,p) 



V'/(0-, P) : 



Wi{crp;a) p-1 p 

with the prime denoting differentiation with respect to the first slot of Wi{z;a). 
Next, we rearrange terms and introduce some new symbols, thereby arriving at 



^^""^ ^G^ + ^^ + l^V N(p)V.(a,p) 



Wl[ap-a) 

ap 



Wi{ap;cT) 

where in terms of the function F{p) appearing in the line element 114(1 we have 
introduced the temporal lapse function N{p) = F-^/^{p) and the radial lapse function 
M(p) = F~^/'^{p). The metrical function N describes the proper-time separation 
between neighboring three-dimensional level-r hypersurfaces, whereas, in a given 
such three-surface, M describes the proper-radial separation between neighboring 



^'^Our notation here is not strict, since of course the RHS has no cr dependence. A strict notation 
would replace each a on the RHS with a blank □ or dot • to indicate that this dependence is lost 
in the integration. 
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p = 1 


P = Pmin 


P — Pmax 


P = P6 



Figure 2. Initial wave-packet configuration in the fre- 
quency DOMAIN. The undulations represent the initial data 
which is compactly supported on [pmin, Pmax]- The domain So = 
(pmax, oo) lies to the right of the data. 



concentric two-spheres |27| . Upon inverse Laplace transformation, the last equation 
becomes 



0-/OS 



Wl{apB;cr) 



Wi{apB]'y) 



with * here indicating Laplace convolution (defined just below). On this equation 
we remark that the direction N^^d/dr + M~^d/dp is null and outgoing, whence 
the derivative of the field appearing on the lhs is along a characteristic. The last 
equation holds in particular at ps, and as our robc wc adopt the following: 

' 1 dijji 1 din in 



(62) 



N Ot 



dp Up 



P=PB 



Pb^N(pb) / uji{t - T';pB)iJi{T',pB)dT' , 



where we have introduced the time-domain radiation kernel (tdrk) 

Wl{apB;cr) 



(63) 



^^i{t;pb) = Cr 



O-PB 



ir). 



Wi{apB;cr) 

We refer to that appearing within the square brackets on the rhs as the frequency- 
domain radiation kernel (fdrk), and we also denote it by 

Wl{apB;cF) 



(64) 



'^;(ct; pb) = crpB 



Wi{apB]o) 
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We assume that the fork i1!i{(t;pb) has the appropriate cr-decay necessary for a 
wen-defined uji{t]pb) = C^'^[uJi{a] pb)]{t). Both uji{t;pb) and u);((t;/9b) do of 
course depend on the values of I and pB (and on the choice of spin j), but to avoid 
clutter we will sometimes suppress this dependence and write simply w(r) and u>{a). 
Finally, we note that the robc can be written simply as 

= Pb'N(ps) I UJl{T-T'-pB)-^l{T\pB)dT' 

in terms of the field '^i = pipi appearing in (|22|) . 

1.4.2. Representation of the kernel. In SECTION |31 we undertake a fairly thorough 
numerical investigation of the analytic behavior of both Wi{apB, cr) and uJiicr] pb) 
as functions of the complex variable a. As a result of our investigation, we shall 
make the following conjectures regarding the fork (1!i{(t;pb)- First, for / fixed 
6ji{a] Pb) is analytic on C\(— oo,0], save for Ni = Ni{pb) € Z>o simple poles with 
locations {ai_n = o'i,nipB) ■ n = 1, • • • , Ni} lying in the lefthalf a-plane. Second, 
Ldi{a] Pb) is bounded in a neighborhood of the origin ct = 0. Third, Rea);((T;ps) 
is continuous and lmuji(a; pb) jumps by a sign across the branch cut along the 
negative Recr axis. The integer Ni{pb) is constant over sizable regions of the pB 
parameter space. However, the pole locations ai^n do vary smoothly with respect 
to changes of ps, apparently subject to 

oo 

(66) 0-l,n{PB) ^ ^ CTi.n.fcPs'' , 

fc=l 

where the cFi.n.k are constants. This series is perhaps only summable in the sense 
of an asymptotic expansion, and we have only numerically observed the first two 
terms. 

Let us now define the nth pole strength and a cut profile respectively via the 
formulae 

(67) aiAPB) = -PB(y'un{PB) , /;(x; Pb) = Imaj;(xe"'; pb) , 

with X ^ and the prime here standing for d/dpB differentiation. Like the pole 
locations, both of these objects also vary with respect to changes of pB as indicated. 
As can be inferred from the third conjecture of the last paragraph, it is the case 
that Im(2);(xe~"^; Pb) = — fiix'i Pb)- To give a concrete example, we choose j = 0, 
I = 2, and pB = 15, in which case we have numerically found that N2{lh) = 2, 
CT2,n(15) ~ -0.0969 ±10.0612, and a2,n(15) ~ -0.0936 ± 10.0647, for n = 1(+) and 
2(— ). For these parameter values the corresponding cut profile is shown in FiG. |31 
The plot is typical in the sense that for all I and pb considered here, /i(x; Pb) decays 
sharply in the x ^ 0^ a-^d X ~* limits (except for I = where the decay in the 
X ^ 0"*" limit is not as sharp). However, the shape of the profile can be qualitatively 
different for other parameter values. Moreover, for certain exceptional values of the 
parameters, the profile can even blow up at a particular x point, in which case 
numerical evidence suggests that the integral in H68|l is defined in the sense of a 
Cauchy Principal Value. We discuss all of these issues in Section 13.11 



(65) 



"1 d^i 
N dr 



1 d^i 
M dp 
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0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

X 

Figure 3. Typical cut profile. For this plot / = 2, = 15, 
and .7 = 0. 



In terms of the pole locations and strengths and the cut profile, we claim that 
the FORK has the following representation (suppressing pB dependence for now): 

68 uji{a) ^ y / — ■ — dx, 

^1 cr - <^i,n TT Jo cr + X 

for all (T not equal to a ai^n and not lying on (— cx),0]. Despite the right closed 
bracket here, we shall also evaluate this representation at cr = 0. Given the de- 
scribed structure of the cr~function a)i(cr), the derivation of such a representation 
amounts to a simple exercise involving the Residue Theorem and a "keyhole" con- 
tour. Although we will often describe the second term on the RHS of H68|) as 
corresponding to a continuous set of poles, these are not really poles in the sense of 
complex analysis (which are properly isolated singularities) . Formally, we compute 
the inverse Laplace transform of H68|l . with result 

(69) uji{t) = V ai,n exp(cr;.„T) / fiix) exp(-xr)dx . 

n=l 

Evidently then, direct numerical construction of the fork would amount to numer- 
ical computation of the pole locations and strengths and the cut profile. For I < 10 
we consider such a direct construction in Section 13.21 although we show in Sec- 
TiON rOl how this brute-force approach may be bypassed insofar as implementation 
of ROBC is concerned. 
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1.4.3. Comparison with the Green's function method. Any temptation to identify 
the pole locations (Ji,n{pB) in the representation (|68|l with so-called quasinornial 
modes j68| should be resisted. For a given I there are an infinite number of quasi- 
normal modes j68j . fixed numerical values intrinsic to the blackhole geometry and 
certainly insensitive to any particular choice of outer boundary radius pB- How- 
ever, for > 1 the poles now under examination, that is the zeros in a of the 
Heun-type function Wi{crpB',<y), are finite in number, and they do depend on pB- 
Moreover, the boundary value problem associated with these Heun zeros is different 
than the usual one associated with quasinormal modes. This usual boundary value 
problem was considered in the pioneering work |38) of Leaver, and more recently in 
a careful study by Andersson j39) . The goal of both authors was to examine a given 
multipole field tl'i{T, p) = p~^4';(r, p) in terms of a Green's function representation 
involving initial data. Andersson refers to this as the initial value problem for the 
scalar field (or, more generally, for electromagnetic or gravitational perturbations), 
although when we mentioned that name in the first paragraph of the introduc- 
tion we did not have this Green's function approach in mind. In Eq. (6) of |39j . 
Andersson expresses the scalar field as 

(70) «'/(t,p,) = J G,(p*,p:;T)a,*,(0,p.)dp:+ J drGiip,,p',;T)-9iiO,pM- 

In this equation we view the field '5; introduced in H44(l as depending on (as 
Andersson does), and we have also slightly modified Andersson's notations to suit 
our own. The appropriate limits of integration in (|7n|l are discussed in ISH] . This 
problem perhaps resembles our own; however, as we now demonstrate, it is different 
both in concept and detail. 

Both Leaver and Andersson considered the (here Laplace) transform of the 
Green's function in (|7()|l . a frequency-domain Green's function G; (p* , p'^ ; cr) as- 
sociated with the following boundary value problem. The solution is pure ingo- 
ing at the horizon [vfj" (cr, p* ) ^ exp((Tp*) as p* —>■ — oo] and outgoing at infinity 
[^^^((T, p*) ~ exp(— CTp*) as p* oo]. In fact, we briefly considered Gi{pt,, p'^,; a) in 
and around H57|l . although we shall make no further use of it in this article or the 
follow-up article. Leaver and Andersson's approach was essentially to examine the 
value ^'/(t, p*), as expressed by (fTD)), via a careful analysis of G; (p* , p'^ ; cr) . (Of no 
concern here, Leaver further considered more general driving source terms beyond 
just the initial data.) When Gi{p^, p'^; a) is considered as an analytic function of 
complex a and continued into the lefthalf plane, its pole locations are the quasinor- 
mal modes and there is also an associated branch cut along the negative Recr axis 
|38[ I39L I65| . These complex analytic features play a prominent role in describing 
the physical behavior of the field (see [HHl OHj and references therein) . 

Our key representation (|68(l stems from continuation into the lefthalf cr-plane 
of the FORK Cji{(J] pb), the expression (|64() involving the logarithmic derivative of 
Wi{ap]a) — exp((Tp*)^j^(cr, p). Here we again view as depending on p, rather 
than p, as in the last paragraph. We stress that the fork uji{a\pB) is not the 
Green's function Gi{pf,p'^;a) considered by both Leaver and Andersson. Indeed, 
0Ji{cr; pb) is a boundary integral kernel. Moreover, it is built solely with the outgoing 
solution to the homogeneous ODE, whereas construction of the Green's function 
requires two linearly independent homogeneous solutions, and As men- 
tioned, the pole locations associated with Lbi{a; Pb) are not the quasinormal modes. 
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rather the special frequencies, finite in number, for which the outgoing solution 
Wi{<7p;a) also vanishes at pB- Despite the fact that Gi(p*, p'^; cr) and Cji{f7\pB) are 
different integral kernels, wc remark that they share the same qualitative features 
in the Icfthalf cr-planc (each has poles and a branch cut). 

On top of these technical differences between our work and those of Leaver and 
Andcrsson, we point out that our overall goal is very different. As mentioned, their 
goal was to examine the actual value of the field via the representation (|7l)|l based 
on spatial convolution. On the contrary, our boundary kernel uji{t; ps) is associated 
with temporal convolution, with the goal being to impose exact radiation boundary 
conditions at a given outer sphere B. That is to say, our goal is domain reduction 
via the introduction of integral convolution over the history of the boundary 
B. With this distinction in mind, compare our key Eq. (|62f) with Andcrsson's 
key equation, as we have written it in (|70|l . Perhaps the approach of Andcrsson 
and Leaver could also be used to numerically implement exact radiation boundary 
conditions in an alternative way [by setting = ps + log(pB — 1) in H7()() ]. but 
they did not address this question per se. Moreover, such an approach would 
necessarily relate robc to the details of the data on the initial surface, which would 
seem awkward from a numerical standpoint. ^"'^ Even were such an implementation 
carried out, memory and speed issues would inevitably arise. Besides developing 
exact ROBC via domain reduction, we also intend to provide an efficient and rapid 
implementation of these conditions. 



1.4.4. Approximation of the kernel. Our numerical implementation of ROBC rests 
on approximation of the exact kernel lo[t) (now suppressing I as well as ps) by 
a compressed kernel C(t). We explain this terminology later, but here collect an 
estimate needed to address the relative error associated with such an approximation. 
We start with the Laplace inversion formula, 

(71) ^(^)^2^ /""^(^)e-da, 

where we are assuming that any singularities of i^i^cr) lie in the lefthalf plane Recr < 
0. A change of variables casts the inversion formula into the form 

(72) ^(^)^^|°°^(iy)ei«My, 

thereby introducing the Fourier transform i^{y) — i^iiy) of iP{t). It then follows 
that 

(73) II^IIl.(O.oo) = IklL.(E) = II^IL.(E)- 

The first equality follows subject to the assumption that iP(t) vanishes for r < 0, 
while the second is Parseval's identity. 

Now suppose ^(r), with Fourier transform ^(y), is an approximation to the ker- 
nel w(t). Later we shall have its Laplace transform ^(cr) as a rational function 



We are definitely not critical of the most excellent works of Leaver and Andcrsson. Indeed, 
just as their Green's function technique would seem not the best way to implement ROBC, we do 
not believe we could directly reproduce their results with our boundary kernel technique. 
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P{a)/Q{(T), and then ^(y) = P{iy)/Q{iy). Also introduce the Fourier transform 
uj{y) = ijjiiy) of ijj{t). In terms of these variables we have 

(74) ||C*^-^*^|L,(o,oo) = II^V^-'^^IIl.(r) 

as our basic estimate. Because of this estimate, we focus on finding approximations 
^(cr) to (2'(cr) which have small relative supremum error along the imaginary axis. 

Finally, suppose that we do not quite know Coiiy). Rather, as is the case, we 
must generate Coiiy) itself numerically. Then, instead of the relative error 

^7^^ sun 

{'^) sup^eiK , 

we should consider an expression like 

I'OJ sup^gjjj h sup gijj— — — j- , 

where the final term is an estimate of the supremum relative error in our knowl- 
edge of Lo{\y). For the methods we develop to generate li(iy), this second term is 
negligible with respect to the first one. 

2. Numerical evaluation of the outgoing solution and kernel 

This section describes the handful of numerical methods used in this work. The 
first subsection describes a numerical method for evaluating the outgoing solution 
Wi[apB]<y) at a given complex tr, and this method allows us to numerically study 
the analytic structure of Wi {apB', cr) as a function of Laplace frequency. The lefthalf 
(T-plane is the domain of interest, and a study of Wi{upB',<j) on this domain, 
carried out in Section 13.11 justifies the key representation Hfci8|l . As we indicated 
in Section I1.4.4I our numerical approximations to the fork Cji{a\ pb) are tailored 
to have small relative supremum error along the Imcr axis. Therefore, insofar as 
implementation of ROBC is concerned, we primarily need numerical methods for 
obtaining accurate numerical profiles for Rea)i(iy; pb) and Ima);(ij/; pb) with y G R. 
The second subsection describes two such methods. 

Before describing our numerical methods, we note that Leaver has analytically 
represented a solution to the Regge- Wheeler equation (more generally to the gen- 
eralized spheroidal wave equation) as an infinite series in Coulomb wave functions, 
where the expansion coefficients obey a three-term recursion relation |44| . Such a 
series can alternatively be viewed as a sum of confluent hypergeometric functions. 
One approach towards our goal of numerically evaluating the outgoing solution 
would be to use the appropriate Leaver series. However, beyond the issue of nu- 
merically solving the relevant three-term recursion relation, numerical evaluation 
of Coloumb wave functions (for complex arguments) is already somewhat tricky 
|60| . Here we describe far simpler methods, which are nevertheless extremely ac- 
curate. Although simple, our methods are very accurate only for a limited range 
of frequencies (which happen to be precisely the frequencies we are interested in). 
The Leaver series is valid over the whole frequency plane. Although we have not 
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compared our methods with the Leaver series, we beheve they are better suited for 
our purposes. 

2.1. Numerical evaluation of the outgoing solution. From now on let us 
simply refer to Wi{apB', cr) as a Hcun function and Wi{apB) as a Bcsscl function. 
This is not quite correct since, as discussed in Section 11.3.21 Wi{apB',<7) and 
Wi{apB) respectively differ from Hcun and Bessel functions by transformations on 
the dependent variable; however, this terminology will streamline our presentation. 
We now present numerical methods for computing the complex value Wi{apB', cr)- 
The methods have been designed to successfully compute the similar value Wi {upb), 
formally W/((tpb;0) in our notation. Wi{<jp) solves the ode 

obtained directly from (|46|1 via the substitution z — a p. Section II. 3. 31 noted that 

Wi{(7p) is a polynomial X]L=o '^"('-'^'°) " '^^ degree / in inverse trp, with coefficients 
c„ given in (|52|l . With the exact form of Wi{(Tp) we could in principle compute 
the value Wi{(tpb) dircctly.^^ Nevertheless, if we pretend that the exact form of 
Wi{(jp) is not at our disposal, then the task of numerically computing Wi{(7Pb) 
shares essential features with our ultimate task of computing Wi{apB] c). The task 
of computing Wi(apB) has been an invaluable model, and for ease of presentation 
we mostly focus on it here. 

2.1.1. Numerical integration. Focusing on the p-dependence of the solution, we 
write Wi{crp) = Yl\i=ni^-n'^~^^)P~^ ■ Since we shall not allow ourselves to evaluate 
Wi{crpB) as X)L=o(''"''^~")/^b"' truncate the series after some fixed number l — p 
of terms, assuming that 

i-p 

(78) W,M^^(c„a-")p-" 

n=0 

is at our disposal. Truncation by hand of this already finite series serves as a 
model for the scenario involving Wi{ap; a), where only a divergent formal series, 
such as the one specified by and (|^ . is at our disposal. With our truncated 
series we can still generate an accurate approximation to the value Wi{apoo), so 
long as Poo is large enough. Let us set poo = scale * pB, with scale a large 
number. Evaluation of the truncated sum and its p derivative at poo then generates 
initial data for the ode. Moreover, the generated data is approximate to the exact 
data {Wi{(Tpoc),Wl'{apoo)} giving rise to Wi{(jp). Here the superscript p denotes 
d/dp differentiation, whereas a prime ' would denote differentiation in argument. 
We stress that our approximation to the exact data can be rendered arbitrarily 
accurate by choosing poo large enough. Finally, we numerically integrate (|77|l in 
p from Poo all the way down to ps, thereby computing a candidate for the value 
Wi{apB)- 

As it stands, the description in the last paragraph is an outline for a stable 
numerical method, provided Recr > 0. However, for the case Rea < of interest 
the described method is not stable. To see why, consider the I = 1 outgoing 
solution Wi{ap) = 1 + (crp)^^ to (|77|1 . As a second linearly independent solution to 

^^Due to the growth 1541 of the Bcsscl coefficients, such direct computation is plagued by 
increasing loss of accuracy as I grows. 
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Figure 4. Two-component path in 2;-plane. The figure de- 
picts the ray-and-arc path described in the text. Here we have 
Pb ~ 20, a = — 1 -I- i, and scale = 4, so that zb = —20 -I- i20 
(marked by an o) and Zoo = 4 * 20 * V2 ~ 113 (marked by an x). 
Typically scale will be much larger, but the value here makes for 
a good figure. 



the I = 1 ODE, take the ingoing solution Zi{ap) = exp(2o'p)[l — [ap)~^]. Further, 
suppose that initial data for the ODE is obtained from a truncated sum as described 
above, with {1,0} in place of exact data {1 + (crpoo)"^, — ''■(cpoo)"^}- The initial 
data {1,0} corresponds to a linear combination aWi{ap) + bZi{ap) with a ~ 1 
and 6 such that hZi{ap^) ~ 0. To fix some realistic numbers, let pB = 20, 
scale = 250 so Poo = 5000, and a = —0.05. Then we compute a ~ 1.0040 and 
bZi{apoo) ^ 8.0320 x 10"^ where b ~ 1.1229 x lO^^^ The exact value we wish 
to calculate is Wi{—1) = 0. However, with the chosen initial data, even an exact 
integration of ^ from poo = 5000 to pB = 20 yields the value 3.0393 x 10^". 
Since error in the initial conditions is exponentially enhanced, the second solution 
Zi(— 0.05p) becomes dominant as p is decreased. 

2.1.2. Two-component path integration. The simple discussion at hand suggests 
that we should complexify the variable p, rotating poo off the real axis by an angle 
large enough to ensure that the product apoa lies in the righthalf plane. Then 
integration along a ray in the complex p-plane from poo towards the complex point 
exp(i0)pB (with ps still real here) would exponentially suppress error in the initial 
conditions. At the end of such a ray integration, a second integration over an arc of 
6 radians would be needed to undo the phase of exp(i0)pB. We effect such a rotation 
of the p coordinate as follows. We choose to work with the variable z = crp, the 
solution Wi{z), and the truncated series X]L~ifo ^n^:"". Our integration will now be 
carried out in the complex z-plane rather than the p~plane, although the strategy 
is essentially the same. We define Zqo to be a large real number scale * |crps|, 
and obtain initial data approximate to {Wi{zao),Wl{zao)} via evaluation of the 
truncated series and its z derivative at z^o- Even for large / we have typically 
chosen I — p = 5 terms to define the truncated series. Then to compute Wi{apB), 
we must numerically integrate the ODE (|46|l from z^o to zb = crpB along some path 
in the complex z-plane. A possible two-component path is shown in FiG.0] It is 
composed of a straight ray followed by a circular arc, with the terminal point of the 
ray being the real z-point \apB\- The arc subtends an angle equal to the argument 
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Figure 5. Mirror image paths. The figure depicts tlie two 
globally different paths connecting z^o to z-points above and below 
the negative Rez axis. Use of such different paths puts a branch 
cut on the negative Rez axis for Heun functions and for Bessel 
functions not of half-integer order. 



of cr. If a happens to lie in the third quadrant, then the relevant two-component 
path looks like the one in FiG. 01 except reflected across the Rez axis. 

Evaluation of the Heun function Wi{apB', cr) features numerical integration of the 
ODE H44|l from Zqo to zb = crpB along the same two-component path. Although 
z of course changes along the integration path, the a in Wi{z]ij) remains fixed 
throughout the integration. We are then integrating a different ode for each value 
of cr. Since these Heun functions are not bispectral (see the discussion in Section 
[T23, there would seem no way around such a cumbersome approach. Were we 
only interested in Wi{apB), and not Wi{(jpB', cr) as well, such an approach would be 
unnecessary (for then we could integrate with respect to frequency a). In essence 
our two-component path method for evaluation of either Wi{apB) or Wi{apB\o) 
is an integration with respect to radius rather than frequency. Indeed, even for the 
Bessel case, we connect each zb to the point Zoo by its own integration path, and 
during the integration do not record values for Wj(z) along the path. Recording 
such values throughout the integration would be a more efficient way of mapping 
out the cr dependence of Wi{<tpb)- 

Let us note two key features of two-component paths. First, for any choice of 
zb ^ the associated path avoids the origin where the function Wi{z) is singular. 
Second, considering two terminal points, one Zg just above and the other Zg just 
below the negative real axis, we note that the respective two-component paths 
connecting them to z^a are mirror images, as depicted in FiG. \S\ Therefore, the 
path leading from Zoo to Zg is globally different than the path leading from Zoo to 
Zg, this being true despite the fact that z^ and z^ may lie arbitrarily close to each 
other in the z-plane. For I e Z>o the functions Wi{z) are clearly analytic on the 
punctured z~plane, so that Wi{zg) = Wi{zg) in the limit that these points meet 
on the negative Rez axis. However, for Heun functions we shall find Wi{zg;(T^) ^ 
Wi(z^;a'^), and in turn Wi{(j^ pB]cr^) ^ Wi{a'° pB]cr'°), for corresponding = 
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Figure 6. Contour lines of log |M^i(15CT)| on ct-plane for 
I = 70 AND 74. On the lhs wc plot log |W7o(15(7)|. The method 
used to generate this plot is the one based on two-component 
paths, and the integration scheme along each path component is 
fifth-order Runga-Kutta-Fehlberg . Relevant parameters here 
are / = 50, J = 50, TV = 49152, M = 49152, / = 70, ps = 15, 
scale — 500, p — 65, and k — 1. I and J respectively specify the 
vertical and horizontal discretization of the tr-plane. N and M are 
respectively the number of integration steps along the ray and arc. 
Other parameters are described in the text. On the RHS / — 74 
rather than 70 and p = 69 rather than 65 (the initial condition is 
still determined hy 5 = I — p terms). All other parameters are the 
same as for the LHS plot. 



^b/pb ^nd <T^ = Zg/pB- Therefore, the negative Recr axis is a branch cut for 
Wi{apB', c) as a function of a. Our path choices for connecting points in the second 
and third quadrants to Zqo have been made precisely to put this branch cut on the 
negative Recr axis. 

As we demonstrate below, the described two-component path method is quite 
accurate for low I. However, for large / and some values of a there is a consid- 
erable loss of precision associated with evaluating Wi{(JPb) by this method (this 
is true no matter what integration scheme is used along the path components). 
Therefore, we shortly introduce a more accurate method based on one-component 
paths. Before turning to the improved method, let us first heuristically describe the 
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Figure 7. Value of log^Q |W74(z)| along two-component 
PATH FROM Zoo TO zr = 15cr FOR a = -3.29128 + iO. 05785. The 
point -3.29128 + iO. 05785 is close to a zero of Wjiilbcr) and lies 
in the region of the cr-plane shown in Fig. The parameters here 
are the same as those listed in the caption of Fig. |H1 Along the 
horizontal axis we have the N + M = 98304 integration steps. The 
function is of order unity at Zoo and close to zero at the terminal 
point zb- However, note that the modulus of Wr4{z) gets larger 
than lO^*' during the integration. 



trouble the two-component method can run into for large I. In FiG. Elwe graph- 
ically demonstrate the breakdown in the method which occurs (for the specified 
parameter values) when / gets beyond 70. The relevant task under consideration 
is to obtain Wi{apB) in a region around those zeros of Wi{apB) which have large 
negative real parts. On the LHS we plot log |W7o(15(t)|, using the logarithm to 
distribute contour lines more evenly. For the portion of the cr-plane shown only 
two of seventy zero locations are evident. Note the onset of degradation in the 
numerical solution. On the RHS we plot log |VF74(15tT)|, and in the plot two of 
seventy-four zero locations are somewhat evident, despite significant degradation. 
This degradation stems from the following phenomenon. Although they do avoid 
the origin, two-component integration paths, especially those which terminate near 
a zero with large negative real part, tend to pass through a region near the origin 
where the solution is quite large. The phenomenon becomes more pronounced as 
I grows. Two-component paths connect Zoo (where the solution is of order unity) 
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Figure 8. Asymptotic curve C. The large curve C is the 
one described in the text and near which the scaled zeros of 
the MacDonald function lie. It has parametric form z(A) = 
— — AtanhAi iV A coth A — A^ for A in the domain [0, Ao] with 
Ao — 1.1997 such that tanhAo ~ I/Aq. The small curve is a circle 
tangent to the C point (xq, 0), where = — \/Aq — 1 ~ —0.6627. 

to zb (which might be at or near a zero of the solution in question), and at each 
of these points the solution is in some sense small. Therefore, loss of accuracy is 
an issue if the connecting path indeed passes through a large-solution region. We 
document an instance of this situation in FiG. [7| 

2.1.3. One- component path integration. We now describe an alternative class of 
integration paths tailored to mitigate the problem of passing through regions where 
the solution is large. Members of this alternative class are one-component paths, 
and this new class yields an improved version of the integration method based the 
two-component paths. As the new method will be more accurate, we will use it to 
quantify the accuracy of the two-component method. 

The one-component paths of interest are essentially dilations of a certain curve 
C depicted in FiG. |S1 A parametric description of C in terms of transcendental 
functions is given in the figure caption. The curve C is intimately related to the 
zeros of Wi{z), also the zeros of the MacDonald function 7^^+1/2(2). As a degree-Z 
polynomial in inverse z, the function Wi{z) has I zeros. Let n € Z>o run from 1 
to I (with ri = if / = 0) and fc/_„ denote the zeros of Wi{z). It is known that the 
scaled zeros {I + l/2)^^fc;_„ He arbitrarily close to C as / becomes large (see results 
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Figure 9. One-component path. To generate the depicted 
curve we have set R ~ 38.7479 and chosen 77 from the range 
0.5162 ^T] < 1.8650, ensuring that the terminal point zb = 20+i20 
and the initial point Zqo ~ 113 + 172 (comparable with the anal- 
ogous point in Fig. Q . For one-component paths the meaning of 
the scale variable is a little different. For such paths the real 
component of Zqo is set by i? * scale. 

listed or summarized in Refs. [611 118L 1631 12 Ip. A simple numerical experiment 
performed in Section 13 . 1 . II confirms this assertion even for small /. Therefore, for 
a given /, dilation of C by Z + 1/2 yields a curve on which the solution Wi{z) tends to 
remain small. Our one-component integration paths are quartic approximations to 
(dilations of) C, and an example is depicted in FiG.El The approximation is given 
parametrically by R{g{r]), rj), where R is fixed and g{r]) ~ arf + hrf' + c is a quartic 
polynomial such that upon multiplication by R the C points (0,±1), (a;i,±2/i), 
and (xo,0) all lie on the parametric approximation. We have xi = — ■\/2/(e^ + 1), 
Ui = ^2/ (e^ — 1), and xo = —\/Xq — I with Aq — 1.1997 obeying tanh(Ao) = I/Aq. 
From the parametric description of C given in the caption of FiG.|Hl one may verify 
that each of these points indeed lies on C. Numerically then a ~ 0.1534, b ~ 0.5093, 
and c ~ -0.6627. 

We repeat the graphical investigation described and carried out at the end of 
Section 12.1.21 but now with the one-component method. The relevant contour 
plots of log 1 1^70(150-) I and log |M^74(15(t)| are shown in FiG. El Comparing this 
set of plots with the corresponding set in FiG. El we see significantly less degra- 
dation in former set. Moreover, we again investigate the size of the solution along 
an integration path in FiG. These figures and their captions argue that the 
method based on one-component paths is more accurate than the one based on 
two-component paths, at least insofar as zero-finding is concerned. 

2.1.4. Accuracy of the numerical evaluation. To check the accuracy of our methods 
for evaluating either Wi{apB', c) or Wi^aps) we may compare values obtained inde- 
pendently from one-component path and two-component path integration. For the 
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Figure 10. Contour lines of log |M^i(15CT)| on ct-plane for 
I = 70 AND 74. On the lhs we have plotted log |W7o(15(t)|. The 
method used to generate the plot is the one based on the one- 
component paths, again with Runga-Kutta-Fehlberg integration. 
Relevant parameters here are / — 50, J ~ 50, P = 98304, I = 70, 
Pb = 15, scale = 500, p = 65, and n = 1. P is the number of 
subintervals for the numerical integration. Other parameters arc 
described in the text or the caption for Fig. The plot on the 
RHS is nearly the same, except that now I = 74 and p = 69 (i.e. 
the initial condition still determined by 5 ~ I —p terms). All other 
parameters are the same as for the lhs plot. 



evaluation of Wi{(JPb) there are other checks. First, numerical values for Wi{apB) 
can be cheeked against direct evaluations of X]L=o Cn(cps)~". However, as we have 
seen in (|54|) . for large / the final coefficient c; and the ones just before it quickly 
become too large to faithfully evaluate this exact expression. One can use extended 
precision (say in Mathematica) to get around this problem. Another check, useful 
for large values of / even without extended precision, involves the known continued 
fraction expansion 

(70) W{iz) _ 1(1 + 1) {l-m + 2) 2(2/ -1) 21 

^ ' ^Wi{z) 2{z+l)+ 2{z + 2)+ 2{z + l-l)+ 2{z + l)' 

This formula follows from recurrence relations obeyed by MacDonald functions [2] ■ 
It remains valid for non-integer I; however, in this case the RHS of the equation is an 
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Figure 11. Value of log^Q |W74(z)| along one-component 
PATH FROM Zoo TO zr = 15(7 FOR a ^ -3.29128 + iO. 05785. The 
point —3.29128 -I- iO. 05785 is close to a zero of H/74(15cr) and lies in 
the region of the cr-plane shown in Fig. EH Along the horizontal 
axis we have the P = 98304 integration steps. Note that the 
modulus of M^74(z) never gets so large as it does in Fig.|7| 



infinite continued fraction. Lenz's method may be used to evaluate this continued 
fraction for any I (see the appendix of Ref. jHOI)- Now, both of our integration 
methods also return the derivative WI{zb) in addition to Wi{zb)- To see why, let 
W ~ U + iV (suppressing the argument and I). In order to integrate the second- 
order ode (|46|I for the complex variable W, we switch to a first-order system of 
ODE for the real vector {U* ,U,V,V*). The * denotes differentiation with respect 
to any relevant path parameter. With knowledge of U* and V* and the Cauchy- 
Riemann equations, one can recover W' . Therefore, both the one-component and 
two-component path methods may also be used to evaluate zbWI{zb)/Wi{zb), and 
this value can then be checked against the continued fraction (|79|l evaluated a± zb- 
In this context, notice that the zeros in a of the reciprocal Wi{apB)l{'JPBWl{apB)) 
are also the zeros of Wi{apB): owing to the fact that the zeros of the MacDonald 
function are simple Appealing to the above checks, we find that even the 

inferior numerical method based on two-component paths is quite accurate for I < 
10; and we offer the following concrete investigations to sharpen this statement."'^'^ 



We remark that these accuracy checks test our methods where we need them most, that is 
on those tasks necessary for a numerical construction the kernel via the representation 1681 . 
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Table 3. Zeros of Wio(15cr) computed via two different 
METHODS. In the top table we list five of the ten zeros (the other 
five arc complex conjugates). These have been found using the 
two-component path method in tandem with the secant algorithm. 
In the second table we list the same zeros, although now found 
using the one-component path method with the secant algorithm. 



Accuracy in zero-finding. The function VI^io(15ct) has ten zeros, which come 
in five complex-conjugate pairs. Using the secant algorithm, we compute the five 
zeros with positive imaginary parts via our two independent methods. Note that 
whether the one-component or two-component path method is used, each function 
call in the secant algorithm involves a numerical integration. The results, listed in 
Table 13 indicate that for low I the two-component path method is associated with 
absolute errors equal to or better than 10"^^, at least insofar as zero-finding is 
concerned. We reach the same conclusion upon computing the zeros of the Heun 
function Wio(15(7; tr) via the two methods. This is remarkable in that there is no a 
priori relationship between the asymptotic curve C and the zeros in a of the Heun 
function Wi{apB',cr). However, carrying out the same graphical experiments for 
Heun functions that wc carried out for Bessel functions and documented in FiGS.[7| 
and 111! we again find that the one-component path method is better than the 
two-component path method at keeping the solution small during the integration. 
Numerical experiments described in Section 13. 1 . II further clarify this issue. 

Accuracy in the cut profile. As applied to the Heun case, both the two-component 
and one-component path methods also return Wj'((tpb; c). This can be seen 
via argumentation similar to that given above in the context of the real vec- 
tor {U* ,U,V,V*). Therefore, we have two independent methods for calculating 
apBWI{c7pB]<7)/Wi{apB',cr), where a may be chosen pure real and negative (say 
with the convention that all paths approach the negative Retr axis running through 
the second quadrant). That is to say, each of our methods may be used to evaluate 
the cut profile 

(80) fiix;PB) - Im [e'^xPBWlie'^xPB;e'^x)/Wiie'^xPB;e'^x)] ■ 



Computing the same zeros in extended precision with Mathematica, we have checked that 
the one-component path method yields the zeros with absolute errors near 10^^^. 
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Figure 12. Absolute error in the cut profile /2(x;15). 
Here we plot the difference A/2(x; 15) of two numerical computa- 
tions of /2(x; 15), one based on the two-component path method 
and the other on the one-component path method. For the two- 
component method we have N = 131072 = M, while for the one- 
component method P = 262144. Parameters common to both 
computations are k = 1, scale = 1000 and p = —3 so I — p = 5. 
There are 512 x-subintervals in the plot. 

Using each method to obtain its own numerical graph for the profile /2(x;15) 
shown in FIGURE 01 of Section I1.4.2I we then plot the difference of these graphs in 
Fig. El Similar graphs for other values of ^ < 10 indicate that the two-component 
path method evaluates the maximum value of |//(x;15)| with an absolute error 
better than 10^^°. For the following reasons we believe that the one-component 
path method computes this maximum with an even smaller absolute error. The 
essential support of |/z(x;15)| corresponds to a region of the a plane near those 
zeros of Wi{apB',cr) with largest negative real part. Therefore, in connecting z^o 
to a purely real zb = cxp(i7r)xpB on the cut, a one-component path runs all the 
way near (a dilation of) C, indicating that the numerical solution along such a one- 
component path again tends to remain small. Experiments like those documented 
in Figs. 171 and II II confirm this expectation. 

We will mainly use the described integration methods for small I < 10. How- 
ever, we note that via comparison with both the continued fraction expression 
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and extended precision calculations in Mathematica, we believe that our one- 
component path method maintains single precision accuracy up to about ^ = 50, at 
least insofar as zero-finding is concerned. Finally, wc mention that wc have carried 
out all integration using the Runga-Kutta-Fchlburg scheme with fixed step-size 
along individual path components. In light of the sufficient accuracy noted here 
and the next subsection, we have not found it necessary to introduce any sort of 
adaptive integration. Furthermore, wc have not found the local truncation error 
estimate (stemming from comparison between the fourth and fifth-order integra- 
tion schemes) provided by Runga-Kutta-Fehlburg to be a useful diagnostic for our 
purposes. Relying on our own accuracy checks, we have simply used the straight 
explicit fifth-order scheme. 

2.2. Numerical evaluation of the radiation kernel. The numerical methods 
discussed in the last subsection work well for values of I < 10, and via H68|l will 
allow us to directly construct sufficiently accurate sum-of-poles representations of 
the FORK. Moreover, even for moderately large I these methods prove useful in 
qualitative investigations of the outgoing solution's analytic structure. However, 
when it comes to building an accurate sum-of-pole representation of the fork for 
high I, the described methods lack the necessary accuracy. 

In this subsection we describe different methods for direct evaluation of the 
radiation kernel itself along the Imcr axis, ones sufficiently accurate even for high 
I. Given accurate profiles for the real and imaginary parts of the radiation kernel 
along this axis, we may then extract an accurate sum-of-poles representation via 
a method described in Section 13.31 and due to Ref. ^Hl- Here we described two 
methods for evaluating loi{(j;pb) when a is pure imaginary, one accurate so long 
as |ct| ^ and the other so long as 7^ |cr| < 1. There is some interval of 
overlap on the Imcr axis on which both methods arc accurate and may be compared. 
We warn the reader that we also use the notation u!i{a;pB) for the Bessel fork 
<ypBWl{apB)/Wi{<7pB)- In order to avoid the confusion which might arise from 
this dual meaning of the symbol c<)i(cr; Pb), hi this subsection we sometimes adopt 
the following notation. For the product of z with the Hcun logarithmic derivative 
we may use 

(81) Wiiz;a) = zWliz;a)/Wiiz;a), 
while for the corresponding Bessel object we may use 

(82) wi{z)^zW!iz)/Wiiz). 

Formally wi{z) ~ wi{z;Q), in parallel with the conventions of Section 11.3.21 
For the Heun case uJi{a] pb) = wi{apB',cr), while for the Bessel case il>i{a;pB) = 
wi{apB)- 

2.2.1. Evaluation of the kernel for large imaginary frequencies. We turn first to 
the evaluation of il!i{a; pb) for tr G iR and \a\ ^ 0. For the remainder of this 
subsection a ~ iy for real y. In this scenario we find it useful to again work with 
the complex variable z = ap. As before, for a given a the terminal evaluation 
point will be denoted by zb = o'Pb, and it lies on the hnz axis. Consider two 
positive real numbers scalei > scale2 and associated z-points zi — scalei + zb 
and Z2 = scale2 + zb- The point zi is analogous to the point Zoo introduced 
before. Further consider a straight path like the one shown in FiG. El running 
through all of these points. Let us now outline the method for obtaining the value 
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Figure 13. Two-component path in z-plane. The figure 
depicts a straiglit path of the type described in the text. Here we 
have pB ~ 20, cr = i, scalei = 20, and scale2 = 40, so that the 
points zi, Z2 and zb respectively correspond to the marked cross, 
square, and circle. Typically scalei and scale2 will be much 
larger, but the values here make for a good figure. 



^i{zb / Pb\ Pb)^ mostly considering only the model Bessel case to streamline the 
presentation. First, using the truncated series X^n^To '^"^r"i compute initial 
values for the ODE (|46|l . Next, we integrate the ode along the straight path from 
zi to 22 (the first portion of the path in FiG. 113(1 . As Rez > along this path, 
we again have exponential suppression of errors both in the initial conditions and 
due to roundoff. The result of this integration is accurate numerical values for 
Wi{z2) and VF/(z2), from which we can directly build a numerical value for the 
kernel Z2WI {Z2) /Wi{z2) at this intermediary point. The assumption here is that Z2 
is still large enough in modulus to ensure that the solution Wi{z2) is not too large. 
Finally, we integrate the radiation kernel itself along the straight path from Z2 to 
Zb, carrying this out as follows. 

Whether we are working with H44|) or (|46|l . we have an ode of the form 

d^d; d$/ 

(83) —^ + R{z;a)—^ + S{z;a)<^>i = 0, 

az^ dz 

so that both wi{z) and wi{z\ cr) obey a first-order nonlinear ode of the form 

(84) '^-^ + ^+R[z-(j)wi + zS{z-a)^Q. 

dz z z 

We use the same symbol wi as the dependent variable here, since as a first-order 
ODE there is only one linearly independent solution. With the accurate value for 
Wi{z2) at our disposal at the end of the first integration leg, we integrate this last 
ODE from Z2 to the terminal value zb (the second portion of the path in FiG. I13f) 
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in order to obtain the desired complex value wi{zb) = zbWI{zb)/Wi{zb)- For the 
Heun case Eq. H84(l is an ode for wi{z;a) = zWl{z;a)/Wi{z;(j), and it is again 
integrated from Z2 to zb, given an accurate value for wi{z2', u). 

The method just described can be unstable if the terminal point zb lies too close 
to the origin. Indeed, notice that some of the terms in the ode (|84() are singular 
at the origin, showing that finiteness of the kernel derivative at the origin depends 
on exact cancellation of singular terms. Round off error will spoil any such exact 
cancellation in an integration towards a terminal zb equal to or near zero. Later we 
demonstrate that the value a)/(0; pb) of the kernel at the origin can be computed 
in closed form for both the Bessel and Heun cases. This raises the possibility that 
the value wi{zb) — or in the Heun case the value wi{zB]<y) — corresponding to a 
non-zero \zb \ <C 1 might be numerically computed via integration of (|84|l out from 
the origin. However, our numerical experiments suggest that this is not a viable 
approach. Moreover, the following analytical reasoning would also seem to dash 
this possibility. For the Bessel / — 2 case we have the exact expression 

(85) ^^^(^)/^^(^)-^^37T3- 

A linearization stability analysis of the Bessel-case ODE (j84|l about this solution 
indicates that small perturbations of zWl^iz) jWiiz) grow exponentially on paths 
running away from the origin along the Imz axis. 



2.2.2. Evaluation of the kernel for small imaginary frequencies. For the reasons 
just laid down, we use a different method for small non-zero imaginary a. The new 
method employs integration in the complex p-plane rather than the z-plane. We 
introduce new positive real numbers scalei > scale2 > ps, a phase factor exp(i0), 
and the following associated p-points (all in polar form): pi = scalei * exp(i0), 
P2 ~ scale2 * exp(i0), p^ = pB exp(i9), and pb- These points define a three- 
component path in the p-plane such as the one shown in FiG.El Let us now outline 
the new method for computing the value cD/(tT; pb), again mostly considering only 
the model Bessel case to streamline the presentation. First, using the asymptotic 
expansion (|78() . we compute initial values for the ODE (|77|l . Next, we integrate 
(|77|l along the straight ray from pi to p2 (the first portion of the path in FiG. I14|l . 
We choose the angle 9 such that Re{ap) > along this path, ensuring exponential 
suppression of errors both in the initial conditions and due to roundoff. The result 
of this integration is accurate numerical values for Wi{ap2) and Wl'{ap2), from 
which we can directly build a numerical value for wi{ap2) at the intermediary 
point p2. Similar to before, the assumption is that p2 is large enough in modulus to 
ensure that the solution Wi{(Tp2) is not too large. Finally, we integrate wi{ap) itself 
along a two-component ray-and-arc path from p2 to the real point ps, by way of 
an intermediate point p^. Such a remaining two-component path is depicted in 
Fig. 1141 as the final two portions of the curve connecting p2 io pB- This integration 
is carried out as follows. 

Whether we are working with 1)43(1 or ((77|l . we have an ode of the form 



(86) 
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Figure 14. Three-component path in p-plane. The figure 
depicts an example three-component path of the type described 
in the text. Here we have pB = 20, scalei = 60, scale2 = 40, 
and 6 = — 27r/5 so that the points pi, p2, ps and ps respectively 
correspond to 60 exp(— i27r/5) (cross), 40 exp(— i27r/5) (square), 
20 exp(— i27r/5) (diamond), and 20 (circle). The depicted path cor- 
responds to a (T value lying on the positive imaginary axis, so that 
between pi and we have Rc(crp) > 0. Better suppression of 
error would be had for 6 = — 7r/2, in which case the portion of the 
path between pi and p3 would lie on the negative imaginary axis. 
However, the final integration would then be over a longer arc, and 
on this final arc we do not expect to have error suppression. There 
seems to be some trade-off here, which is why we have kept as a 
parameter. In any case, typically scalei and scale2 will be much 
larger, but the values here make for a good figure. 

whence both wi{(Tp; a) and wi{ap) obey a first-order nonlinear ode of the form 

(87) ^ _ !^ + !i + 7^(p; a)wi + pS{p- a) = . 

dp p p 

To reach this equation, we have used, for example, wi{pa) ~ pWl'{pa)/Wi{pa). 
Given the initial value for wi{ap2) obtained from the first leg of the integration 
in the last paragraph, we first integrate (|87|) along the straight ray from p2 to ps 
which has the same modulus as the terminal point pB- The final leg, a rotation 
back to the Rep axis, is an integration of (|87|l along an arc from p3 to ps. 



40 




-15 -10 -5 5 10 15 

y 




-15 -10 -5 5 10 15 

y 



Figure 15. Bessel fork t2)64(iy;15) = W64(15iy). Here we 
plot the functions it64(15ij/) = Re?i'64(15ij/) and Ve4{15iy) = 
Ini?i;64(15i2/), with the y axis split into 512 subintervals. For \y\ > 
break = 1 we have evaluated 0)51 (iy; 15) using two-component 
integration with the following parameter values: N = 131072, 
M = 131072, scalei = 1000, and scales = 100. N and M 
are respectively the number of integration steps taken along the 
first and second components of the path. For |y| < break wc have 
evaluated 0)51 (iy; 15) using three-component integration with the 
parameter values N = 131072, M = 131072, P = 2048, 9 = 7r/4, 
scalei = 1000, and scales = 100. N, M, and P are respectively 
the number of integration steps taken along the first, second, and 
third components of the path. For both integration methods k = 1 
and p = 59. Typically, wc have chosen break smaller, but now have 
break = 1 to demonstrate the three-component method. 

2.2.3. Value of the kernel at the origin. For the Bessel case the origin value (0; pb) 
of the radiation kernel is the limit Imia-^o 'Wi{apB), while for the Heun case the 
value cD/(0;pb) is the limit liiaa^oWi{crpB',(y)- Whether considering the Bessel or 
Heun case, we may derive an exact expression for the value uji{0;pb)- Turn first 
to the Bessel case, where Wi{ap) = X]L=o c-n{.cr p)~'^ is of course singular at tr = 0. 
However, with this exact expression it is easy to check that 



(88) 



lim wi{(jpb) = -I 

IT— »0 
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Figure 16. Error in Bessel fork a)64(iy;15) = weiilBiy). 
Here we plot the absolute error | Ai(;64(15iy)| as well as the rel- 
ative error |Aw64(15iy)|/|w64(15i?/)|. These errors have been com- 
puted agamst the "exact" W64(15iy) generated with the continued 
fraction expression (|79|l . Parameters are the same as in FiG. 1151 



We stress that this calculation of (l!i{Q; ps) makes use of the exact form of the 
outgoing solution, which is not at our disposal in the Heun case. A separate recipe 
for getting this value, one without appeal to the exact form of Wi{apB), goes as 
follows. Set cr = in H77|l. thereby reaching an ODE 



(89) 



dp^ p2 



with solutions and p We now use [pdplogp ']|p=pB as the origin value 
uji{0;Pb), again finding -I. 

Let us turn to the Heun case and follow this recipe for getting the value uji(0; ps)- 
We set (T = in obtaining the following ode: 



(90) 



d^l., 


■ 1 1 ■ 


d$j 


'k K + l{l + l) 


dp2 + 


p p-l 


dp + 


_p2 p{p-l) 



= 0, 



42 



X 10" 




1000 



200 400 600 800 1000 



Figure 17. Error in Bessel fork L2;64(iy; 15). These are es- 
sentially the same plots as those in FiG.El save that here we have 
a larger y-interval. 



Both solutions to this equation may be expressed in terms of infinite scries in inverse 
p. The one corresponding to above has the form 



(l+n) 



(91) 

where oq = 1 and 

^^^^ " \l + n+l){l + n + 2)- l{l + 1) ■ 

The series is positive and absolutely convergent for all p > 1. Wc then have 



(/ + n){l + n + 2) + K 



(93) 



wz(0;ps) = - ^{l + n)anPB'' 



n=0 



ri=0 



as our concrete expression for the value in question. Notice that this value ap- 
proaches ~l in the pB ^ oo limit as expected. 

2.2.4. Accuracy of the numerical evaluation. FiG. llSl dcpicts the real part W64(15ij/) 
and the imaginary part ti64(15iy) of the Bcsscl fdrk ct>64(iy; 15) — i064(15iy) along 
the Imcr axis for I = 64 and pB ~ 15. We have generated these plots using the 
methods described in this subsection, and have listed other parameters set while 
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Figure 18. Heun fork a)64(iy; 15) — WG4{15iy- iy). Here we plot 
the functions Mg4(15ij/; ij/) = IlcwQji{15iy; iy) and VQ4{lbiy;iy) = 
Ini'u;64(15i?/; iy). All parameters in these plots match those listed 
in Fig. El depicting the Bessel fdrk. 

obtaining them in the figure caption. To examine the accuracy of these numerical 
profiles, we may compare them with corresponding profiles obtained via the contin- 
ued fraction expansion H79|l . We consider the profiles stemming from the continued 
fraction expansion as the "exact" ones. With the two sets of profiles, one may 
compute corresponding absolute and relative error measures. We plot these errors 
in Fig. 1161 and FiG.EI(the second being a pull-back of the first). From these fig- 
ures we conclude that our numerical methods evaluate WQ4{15iy) with an absolute 
supremum error less than 10~^^ and a relative supremum error less than 10~^^, 
at least for \y\ < 1000. For the Bessel case at hand we have found comparable 
error bounds associated with all other values of I E {10, 11, . . . ,64}, although we 
note that the corresponding y-interval needs to shrink by as much as an order of 
magnitude to maintain these bounds for I = 10. 

Fig. 1181 depicts the real part UQ4(15iy;iy) and the imaginary part 7j64(15iy; iy) 
of the Heun fdrk ii;64(iy; 15) = W64(15iy;iy) along the Imcr axis. We have again 
chosen the representative case I = 64 and pB = 15, setting the rest of the parameters 
to the same values used to generate the Bessel profiles depicted in Fig. El Note that 
the two sets of profiles are qualitatively very similar. However, they are different. 
In particular, now the real part has a minimum value of —66.2816976576098 rather 
than —64. For the Heun case at hand we have no analog of the continued fraction 
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Figure 19. Error in Heun fork (2'64(i?/; 15) = we4{15iy; iy) . 
Here we plot the absolute error measure | Ai(;64(15i2/; and the 
relative error measure |Aw64(15i?/; i?/)|/|'u;64(15i2/; iy)! described in 
the text. We have y £ [0.5, 8] for both. All other parameters set in 
generating these plots are the same as those listed in the caption 
of FiG.ini 

expansion with which to check the accuracy of the profiles. Nevertheless, at least 
for y values of order unity, we can perform an accuracy check by comparing the 
two-component and three-component path methods for evaluating the kernel. Such 
a comparison is shown in FiG. 1191 With the two numerically obtained kernels we 
form an absolute error measure |Aw64(15iy; iy)| and also a relative error measure 
|Au)64(15iy; iy)|/|w64(15iy; iy)|, over y € [0.5,8] for both. In the denominator of 
the relative error, we happen to have used the kernel stemming from the three- 
component method. FiG. ^1 displays plots of both error measures. Note that 
poor performance for the three-component method is evident in the right portions 
of the plots. For the three-component method the length of the third and final 
integration path grows with y. Therefore, for large y one expects a corresponding 
loss of precision for the three-component method. 

3. SUM-OF-POLES REPRESENTATION OF THE RADIATION KERNEL 

In this section we focus on both exact and approximate representation of the 
FORK u}i{(t;pb) as a sum of poles. In the first subsection we qualitatively discuss 
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the exact representation (|68|l of the fork as a (continuous and discrete) sum of 
poles, highhghting what we beheve to be its main features. In the second subsec- 
tion wc document our particular numerical construction of the fork as a sum of 
poles, and give an analysis of its numerical error. We stress that from a theoretical 
standpoint we are conjecturing that the Schwarzschild fork — built from the Hcun 
function Wi{z]a) — admits the representation H68(l . although wc do provide com- 
pelling numerical evidence for a representation of this form. This is in contrast to 
the case of the flatspace fork — built from the Bessel function Wi{z) — which we 
theoretically know admits such a representation |18| . Since for us the representa- 
tion (|68() of the Schwarzschild fork is ultimately conjecture, there is more need to 
painstakingly justify it numerically, and we do so in the second subsection. In the 
third subsection we turn to kernel compression, by which we mean approximation 
of the FORK by a proper rational function P{a) / Q{a) which is itself a sum of poles. 
In the first and second subsections we consider only the j = (k = 1) case, but 
also consider the ] = 2 {k — —3) case in subsection three. As yet, wc have not 
examined the j = 1 (k = 0) case corresponding to electromagnetic radiation. 

In this section and in the simulations we describe later, we have almost exclu- 
sively worked with an outer boundary radius pB G [15,25], corresponding to a 
physical outer boundary radius G [30m, 50m]. Therefore, were we considering a 
more general isolated source of gravitational radiation, one with gravitational radius 
2m, then the boundary two-sphere B would be located outside of the strong-field 
region as defined by Thorne |64j . Moreover, for wave simulation on a fixed back- 
ground as we consider here, the location of B corresponds to a metric coefficient 
F{pb) from H14|l in the range 0.933 < F{pb) < 0.96. Whence B lies in a region 
where the Schwarzschild metric is flat up to small correction. The robc described 
in Section 11.41 are not tied to the weak-field region. However, for this region our 
numerical methods for examining/constructing the FDRK are accurate. 

3.1. Qualitative study of pole locations and cut profile. Using the one- 
component path method described in Section 12.1.31 we first turn to the analytic 
structure of Wi{apB] cr) as a function of frequency a in the Icfthalf plane, assuming 
that a zero of this function corresponds to a pole location appearing in (|68|l . Using 
the same method, we then draw some quick observations concerning the cut profile 
fi{x] Pb) hi (|68|) . We mainly focus on the restricted parameter space § determined 
by Pb G [15, 25] and < Z < 10, but also make mention of some remarkable features 
which crop up for other parameter values outside of this space. Our parameter space 
§ has been chosen with the following reasons in mind. First, its pb interval is as 
discussed in the last paragraph. Second, it includes the first few values of /, which 
we want to single out for special attention. Third (and related to the first two), 
it avoids by design the aforementioned remarkable features. We stress that our 
discussion in this first subsection is mostly qualitative and amounts to a collection 
of conjectures without substantial numerical or analytical proof. Although we are 
bypassing a truly thorough study of some interesting phenomena, we do not believe 
these phenomena to be directly relevant for numerical implementation of robc 
(further remarks on this point to follow). 

3.1.1. Zeros of the outgoing solution as a function of a. Recall that in Section 
12.1.31 wc denoted by {ki ^ : n — \, - ■ ■ ,1} the zero set of the MacDonald function 
Kij^ij2{z) which is also the zero set of Wi{z). With this notation the zeros in a of 
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Figure 20. Scaled zeros of MacDonald functions. Here 
we plot scaled zeros (l + l/2)^^fc;.„ for I = 1,2,3,4. The cross is 
the scaled zero of Ki/2{z), the diamonds are the scaled zeros of 
-^3/2 (-z) I the circles are the scaled zeros of K^/2{z), and the stars 
are the scaled zeros of K-^/2{z). To the eye these zeros, correspond- 
ing as they do to small / values, already lie close to the curve C 
shown and described in both the text and the caption of Figure 
IHlin Section f2. 1.31 As / gets large the scaled zeros (l + l/2)"^fc;.„ 
lie closer and closer to C. 

Wi{apB) are then simply the h^n/pB- Let us collect several facts concerning such 
sets, summarizing results derived or listed in Rcfs. [fill 1181 I21| . First, for even 
I these zeros come in complex-conjugate pairs, while for odd I they again come in 
complex-conjugate pairs save for a lone zero which lies on the negative Recr axis. 
Second, the scaled zeros {l+l/2)^^ki^n lie close to the asymptotic curveC introduced 
in Section 12. 1.31 See FiG.[5D|for a graphical demonstration of this claim. Hence, 
for each I one may imagine the zeros distributed in a crescent pattern in the Icfthalf 
(T-planc. As concrete examples, the zeros of W2{15a) are approximately — O.lOOOi 
i0.0577, while those of WailBa) are approximately -0.1226±i0.1170 and -0.1548 + 
iO. Respectively, these zero sets arc marked by crosses in the lefthand plots of 
Figs. El and m 

Zeros for chosen parameter space. Dealing with Wi{apB', cr) as a function of a in 
the Heun scenario, we have denoted the zeros of this function by ai^n = o-i,n{pB)- 
Over the range S of parameters mentioned above, the zeros ai^n behave qualitatively 
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Figure 21. Zeros of W2{l5a;a) and VK2(25CT;cr). The left- 
hand plot shows contour hnes of log |W2(15cr; cr)|, and the right- 
hand one shows contour lines of log IW2 (25(7; cr)|. The logarithm 
mollifies the singularity at the origin, distributing contour lines 
more evenly. Notice that the zero locations are closer to the ori- 
gin in the righthand plot. In the lefthand plot and the right- 
hand plot we have also respectively marked as crosses the zeros 
of the Bessel functions M^2(15cr) = 1 + 3{15a)-^ + 3(15cr)-2 and 
W2(25cr) = 1 + 3(25cr)"i + 3(25cr)"^ Perhaps evident even to 
the eye, the Bessel and Hcun zeros lie closer to each other in the 
righthand plot (corresponding to the larger value oi pb)- 

similar to the zeros ki^n/pB of Wi{apB), save for one key difference associated 
with odd I. Fig. [2] displays contour plots showing zero locations for W2{l5a;a) 
and W^2(25f7;CT). Likewise, FiG. |23 displays contour plots showing zero locations 
for M^3(15cr; cr) and W^3(25cr; cr). These plots exhibit the main features associated 
with the zeros of Wi{apB',(^) over S. To the eye, apparent zero locations in these 
plots nearly match zero locations (marked as crosses) for the corresponding Bessel 
functions. However, as shown in FiG. 1221 and also noted in the figure caption, there 
is a key difference associated with odd I. To appreciate the difference, compare 
the zeros of W3(15ct) with those of 1^3(15(7; a). Notice that the single zero of 
W3(15(t) lying on the negative Rea axis corresponds to two zeros of W3(15cr;cr), 
one lying just above and the other just below the negative Rea axis. This is a 
generic feature belonging to all odd values of I € {1,3,5,7,9} and pB & [15,25] 
considered here. Therefore, for the Heun scenario and the parameter space § at 
hand there are an even number of zeros whether I is even or odd. Moreover, we 
believe that this feature (of a pair of complex-conjugate zeros lying close to the 
real axis and together corresponding to a single Bessel zero) persists as pB gets 
large, as evidenced by FiG. 1231 and its caption. More precisely, if Ni = Ni{pb) 
denotes the number of zeros belonging to Wi{apB', cr), then for each I £ {0, ... , 10} 
we observe that Ni is constant on [15, 25] with Nq = 0, Ni ^ 2 ^ N2, N3 = 4 ^ N4, 
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Figure 22. Zeros of W3(15cr;cr) and W3{25<t;<t). Thclcfthand 
plot shows contour lines of log 11^3(15(7; a)\, and the righthand one 
shows lines of log |VF3(25cr; cr)|. Again, crosses mark the corre- 
sponding Bessel zeros. Notice in the lefthand plot that the single 
real Bessel zero corresponds to a pair of zeros in the Heun case. 
Actually, in the righthand plot there are also two distinct Heun 
zeros, each one near the Bessel zero lying on the real axis, but the 
resolution is almost too low to see them as distinct. 

N5 = 6 = Ne, Nj = 8 = Ns, and iVg = 10 = Niq. Our conjecture is that each Ni 
is also the same constant on [15, 00). 

Asymptotics of zeros. We discuss two asymptotic regimes for zero locations: one 
large at fixed / and the other large I at fixed ps- Turning to the first regime, we 
conjecture that the zeros <Ji^n{pB) approach the Bessel zeros ki^n/pB as ps becomes 
large. That is to say, the first term in the asymptotic expansion H66() is 

(94) ai^ns ki^n- 

FiG.l^is a typical piece of graphical evidence indicating such behavior. Using the 
zero locations shown in this figure, we have confirmed for each n that \<Tin.n{pB) — 
fcio.Ti/psI = 0{p~j^) over [15,25], in parallel with the first two terms in (|^ . 

Now turning to the second asymptotic regime, we remark on the order scaling 
of Heun zeros as I becomes large. We have observed that the scaled zeros (/ + 
l/2)~^(ji^n tend to accumulate on the same fixed curve Cp^ as I becomes large. For 
example, using the 20 scaled zeros a2o,n(15)/20.5 to hint at a candidate C15, one 
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Figure 23. Zeros of W3(100ct; a). Here we plot contour lines of 
log |W3(100(t; (t)|. Note that two zeros of M^3(100(t; cr) lie outside 
of the plot, as we focus on the pair of zeros closest to the real axis. 
The red cross is a zero of W3(100cr). We might gather from this 
plot that the feature associated with odd / and discussed in the 
text persists as gets large. 

finds that — to the eye at least — all scaled zeros for / < 20 lie on this curve. Of 
course Cp^ is of a different shape than the dilated curve C/ps, although choosing 
a still larger fixed ps value yields better agreement between the two curves. That 
is to say, consistent with the first type of asymptotic behavior discussed, wc have 
PsCpg — !• C holding pointwisc as ps oo. Wc are not confident that this order 
scaling is robust for very large /. 

Zero pair creation as pB approaches unity. We have claimed that over the pB 
interval [15, 25] Heun zero behavior is similar to Bessel zero behavior, save for the 
aforementioned curious feature concerning odd I. As mentioned, we believe our 
claim remains true over [15, oo) for all low I here of interest. However, as we now 
argue, the number Ni of zeros for a given / is not conserved as pB 1^- The 
behavior wc have noticed is the following. 



l^The Heun zeros under consideration are analogous to the "flatspacc quasinormal modes" 
discussed in the introduction of IbSI by NoUert and Schmidt. Nevertheless, despite the infinite 
number of zeros as ps — ► l"*", they are apparently not what physicists call "quasinormal modes." 
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Figure 24. Zeros of Wio(crpB; ct) and Wioicrps)- Here wc plot 
zeros of these functions for pB = 15,17,19,21,23,25. The outer- 
most crescent of diamonds are the zeros of Wio{15a; a), while the 
innermost crescent of diamonds are the zeros of Wio{25a; a). The 
outermost crescent of crosses are the zeros of Wio{15a), while the 
innermost crescent of crosses are the zeros of Wio (25(t) . 



For each value / e {0, . . . , 10} there is a critical value pg of pB (less than 15 
of course) for which a new pair of zeros is created on the negative Retr axis, the 
branch cut. Below that there is yet another critical value for which yet another 
new pair of zeros is created on the negative Retr axis, and so on. Therefore, it 
would seem that as a function Ni{pb) is step-like and blows up as ^ 1+. Let 
us remark on the nature of the zero lying on the branch cut for a critical value 
Pg. As /9b is increased past , say, two zeros appear to collide on the branch cut. 
However, they do not merge into a double zero, rather they pass "over and under 
each other," with each zero remaining on its own analytic neighborhood continued 
across the branch cut. 

In Figs. I?5l and 1261 we document the creation of the first and second new zero 
pairs for W2{apB', cr) associated with decreasing ps from 4 down to 1.5. In Table^l 
we have listed approximations to the critical values for I G {0, . . . , 10}. Using 
<y = X cxp(i7r) along the negative Reu axis, the table also lists approximate values 
for the location x'^^ of each zero-pair creation. We have obtained these numbers 
using the method discussed below in reference to FiG. 123 AH of the approximate 
critical values in the table are well below our pB interval [15,25]; however, for 
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-1.2 -1 -0.8 -0.6 -0.4 -0.2 -1.2 -1 -0.8 -0.6 -0.4 -0.2 

Reo Reo 

Figure 25. Zeros of W2{ctpb;<j)- On the left wc have a plot 
of log |W2(2(t; (t)| and on the right one of log |W2(4cr; a)\, demon- 
strating that a new pair of zeros is created between ps — ^ and 
Pb ^ 2. The initial stages of this creation are perhaps evident in 
the rightmost plot. 



I 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 




1.49 


1.52 


2.52 


2.06 


3.81 


2.67 


5.11 


3.31 


6.44 


3.97 


7.75 




0.63 


1.31 


0.79 


1.49 


0.85 


1.60 


0.89 


1.67 


0.91 


1.71 


0.92 



Table 4. Approximate critical values of ps- Here we list 
rough values corresponding to the creation of the first new zero 
pair for Wi{apB', cr). For example, as pB is lowered from 1.5 to 
1.48, the number A'o of zeros for Wo{(tpb;(j) jumps from to 2. 
We also list approximate values for the location — x'^^ of each zero- 
pair creation, again with ±0.01 error bounds. 



higher I values creation of a zero pair can occur in our interval. For instance, in 
what follows we determine that creation of the first new zero pair for I = 22 occurs 
for p'j^ ~ 15.70 and x'^' ^ 0.96. 

3.1.2. Parameter dependence of the cut profile. Let us now discuss the cut profile 
fiixiPs) appearing in the representation of the FORK. We first remark on 
the behavior of the profile over the chosen parameter range §, and then turn to 
exceptional behavior associated with critical parameter values lying outside §. 

Behavior over the chosen parameter range. For < ^ < 10 and for pB = 15 
and 25 wc plot scaled even profiles in FiG. 1271 and scaled odd profiles in FiG. EHl 
The scaling allows us to view all profiles on the same plot. Notice that the order- 
scaling is different for even and odd cases. As pB is increased towards 25, the other 
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Figure 26. Zeros of 1^2(o'/3b; ct). On the left wc have a plot 
of log |M^2(1.5o-; cr)| and on the right one of log |M^2(1.75(t; (t)|, 
demonstrating that a second new pair of zeros is created between 
PB = 1.75 and pB = 1.5. Asps^l^ more and more zero pairs 
are created. 

endpoint of our interval, all of these profiles retain their shape; however, both their 
maximum value (in absolute value) and the essential window of their support vary. 

Cauchy principal value. The chosen parameter space § has been carefully tailored 
to avoid the exceptional situation where a zero lies on the branch cut. However, 
over our pB interval [15, 25], we shall of course be interested in I values higher than 
10, and on this interval zero pair creation is an issue for such /. A glance at the 
form H8U|) of the cut profile fi{x', Pb) given earlier indicates that a negative real zero 
a = exp(i7r)x of Wi{apB', c) should give rise to a singular cut profile. 

For the aforementioned exceptional case / = 22 and ~ 15.70, we depict the 
profile blow-up in FiG.[52| In the top plot we have the cut profile /22(x; 15.695964), 
where 15.695964 is approximately the critical value of pB corresponding to the 
creation of the first new zero pair for the function W22{<JpB', c). For ps values larger 
than the function has 22 zeros, but as pB is lowered below a new pair of 
zeros appears from the branch cut. The lower plot depicts ReVF22(15.69 5 9 64cr; ct) 
(solid hne) and ImW22(15.69 5 9 64(T; ct) (dotted line) as well as their intersection 
point below the zero line. For pB = 15.695962 this intersection point lies above 
zero, while for pB = 15.695966 it lies below zero. As this intersection point appears 
to move smoothly with varying p^, we conjecture the existence of a zero on the 
branch cut for a critical p^ between pB = 15.695962 and ps = 15.695966 (actually 
we know it lies between ps = 15.695962 and ps = 15.695964). Our guess at the 
value, Pb ~ 15.695964, should be within 2 x 10^^ of the true p^. Furthermore, we 
note that for pB slightly above the critical value the profile /22(x; Pb) is a positive 
peak like one in FiG. [23 but as ps is lowered past pg the profile transitions to a 
negative (and sharper) peak like one in FiG. EHl 
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Figure 27. Even cut profiles scaled by order for pb = 15 
AND 25. Here, for example, /o(x;15)/0.5 is on the far left, and 
/io(x; 15)/10.5 on the far right. 



Despite the blow-up discussed in the last paragraph, we emphasize that along 
the Ima axis, the fdrk uJ22{iy',PB) itself changes smoothly as pB varys across p^. 
Indeed, the pieces Reil'22(iy; Pb) and Ima)22(i2/; Pb) may be computed either via the 
representation H68() or numerically via the methods outlined in Section I2.2I Using 
the latter methods, we observe that both pieces vary smoothly as pB varys across 
the critical value . We therefore offer the following conjecture. Although the cut 
profile /22(x; Pb) is singular at a particular point Xci — 0.96 when pB = p'^b — 15.70, 
the corresponding integral contribution 

(95) _irM^dx 

Jo W + X 

to uj 22 {^y, Pb) varys smoothly as pB varys across pg. This would seem to indi- 
cate that /22(x; PB)/(^y + x) has a distributional interpretation, and we believe its 
integral to be defined in the sense of Cauchy Principal Value. The nearly antisym- 
metrical blow-up in the top plot of FiG. 123 is in accord with this conjecture. 

3.2. Numerical construction of the radiation kernel. We now document our 
numerical construction of the representation (|68|l over the chosen parameter space 
§, also discussing in detail the accuracy of the construction. 
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Figure 28. Odd cut profiles scaled by order for pb = 15 
AND 25. Here, for example, /i(x;15)/1.5 is on the far left, and 
/9(x; 15)/9.5 on the far right. 



3.2.1. Construction for chosen parameter space. To obtain the numerical kernel, 
wo have used the one-component path method described in Section 12. 1.3l both to 
obtain pole locations and strengths as well as the cut profile for certain parameter 
choices. Let us first discuss our treatment of the poles. 

Construction of pole locations and strengths. Let a choice of I G {0, . . . , 10} 
remain fixed throughout this paragraph. We choose nine Chebyshev points {l/p% ■ 
A; = 0, 1, . . . , 8} on the interval [1/25, 1/15]. That is to say, the formula 

(96) 2^1 = 0.7 + 0.04 + (0.7 - 0.04) cos[7r(2A: + l)/(2n + 2)] . 

determines the nine numbers = i/pg. Our choice of nine Chebyshev points 
suffices for our purposes, although we make no claim that nine is the optimal 
number of points. Next, for each k we have used the secant algorithm to find 
the zero set {(yi,n{PB)} of Wi{apg; a). Then, at fixed I and n we interpolate each 
function ai^nips) by an eighth degree Chebyshev polynomial Ti^„{^/pb) in inverse 
Pb, so that this polynomial approximates (Ji^n{pB) on the interval [15,25]. On the 
same interval, we approximate the pole strengths Q;/.ri(ps) by Tl ^{1/ pb)/ Pb, where 
here the prime ' denotes d/d^s differentiation. 

Construction of the cut profile. Given any small positive 77, say ~ 10^^^, we 
assume the existence of corresponding finite integration limits, Xmin (which may or 
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Figure 29. Blow-up in cut profile near a critical value. 
The value 15.695964 is close to the critical value for I = 22. 
See the text for more description. 



(97) „ir-^tePB)dx 



may not be 0) and Xmax: such that the integral 

1 fi{x;p- 

approximates the true value 

(98) _irAtePB)dx 



Jo W + X 

to within relative error rf uniformly for y g R. We stress that this is an assumption, 
although one apparently true for the analogous cut profiles stemming from integral 
order {y = I + 1/2 = n) Bessel functions VF(„_i/2)(o'Pb), as shown in the fourth 
section of ^Sl- In practice we have "eyeballed" the integration window [xmim Xmax], 
for example referring to FlGURE|31of Section I1.4.2I we have chosen [0.0005, 1.125] 
for / = 2 and pB = 15. The correctness of our guess will be confirmed when we later 
examine the accuracy of the kernel. Finally, to obtain the cut contribution to the 
value a)/(iy; pb) for a given y, we discretize the integral (|^ via the Simpson rule. 
For the profile in the aforementioned figure we have used 1024 subintervals. Values 
fiiXj'^PB) belonging to nodes Xj in the corresponding discrete sum are computed 
with the one-component path method. Since for our chosen parameter range the 
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Figure 30. Error in Heun fork d)io(i?/;15) = wio{15iy;iy). 
With 256 subintervals of the y direction, we plot the absolute error 
I Aw;io(15ii/; iy)\ and a relative error |Ai(Jio(15i?;; i?;)|/|wio(15iy; iy)| 
described in the text. 

considered profiles appear to be of definite sign, we expect that these sums are not 
plagued by cancellation error. 

The foregoing construction of the cut contribution is applicable for a fixed value 
of pB; whereas our construction for poles yielded locations and strengths over the 
whole pb interval [15,25]. In order to handle the cut contribution over the whole 
interval, we again introduce the nine Chebyshev points {1/ p% ■ fc = 0, 1, . . . , 8}, and 
— for each Xj integration node — construct an eighth degree polynomial in 1/pb 
which interpolates fi{^5xj / Pb] Pb)- Notice that we are also scaling the integration 
nodes Xj associated with = 15 [and determined by the choice of Xmin a-nd 
Xmax as well as the number of subintervals chosen to evaluate the integral H97|l via 
Simpson's rule]. 

3.2.2. Accuracy of the construction. We check the accuracy of our numerical kernel 
in two ways, and as one result provide compelling numerical evidence that the fork 
indeed admits the sum-of-poles representation 

Value of the kernel at the origin. Building the pole and cut contributions to 
the kernel as described, we may compute a numerical value for (jJi{0; Pb) and check 
it against the accurate series (|?^ . We find that for any choice of I and in § 
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Figure 31. Error in Heun fork (213(12/; 15) = W3(15it/;iy). 
Again with 256 subintervals of y, we plot the absolute error 
|A?i;3(15i?/; i?/)| and a relative error | Aw3(15iy; iy)|/|w3(15iy; iy)| 
described in the text. 



the numerical value for w;(0; Pb) has absolute error less than 10~^^ (in fact on the 
order of 10~^^ or better). We stress that this level of accuracy holds even when 
1/ Pb is not a Chebyshev node, in which case the pole and the cut contributions to 
the kernel are obtained via interpolation. 

Direct check. For \ <l < IQ we now have two independent numerical methods for 
evaluating the Heun fork (l)i{iy\ pb)- The first is evaluation of the numerical kernel 
directly constructed via the representation H68(l as described in this subsection. The 
second is evaluation using path integration as described in Section I2.2I As a final 
and perhaps most convincing accuracy check, we compare these two methods. In 
Fig. 1201 we have such a comparison for a)io(iy; 15) = u;io(15iy; iy). With the two 
numerically obtained kernels we form an absolute error measure | Az/;io(15i?/; i?/)| 
and a relative error measure |At(;io(15ij/; ij/)|/|i/;io(15iy; iy)|, where in forming the 
denominator of the relative error we happen to have used the directly constructed 
kernel. We plot these error measures in the figure. FiG. 1311 depicts similar plots for 
|Aw3(15iy; iy) \ and |Awio(15iy; iy)|/|w3(15iy; iy)\. Over all of §, save for / = cases, 
this check indicates that we have relative and absolute errors better than 10~^^. As 
for ^ = 0, we believe the integration methods of SECTiON l2.2l to be less reliable than 
the directly constructed kernel. Indeed, uJoivy; pb) is quite concentrated around the 
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origin, and very small values of y negate the exponential error suppression built 
into the three-component integration method of Section 12.2.21 In any case, for 
I = the first accuracy check of d)o(0;pB) should be convincing in and of itself. 
Indeed, in the absence of a pole contribution to the kernel, one expects the largest 
error at the y-origin. Moreover, the largest error is indeed concentrated near the 
origin in other small-/ plots such as those shown in FiG. 1311 Via comparison with 
the series (j93|l , we have found that our directly constructed numerical kernel yields 
a value for d)o(0; 15) with an absolute error better than 7.5 x 10^^'^ and a relative 
error better than 2.2 x 10^^^. 



3.3. Approximation of the kernel by rational functions. Throughout this 
subsection we suppress all I and ps dependence (and we continue to suppress k = 
1 — spin dependence as always). 

3.3.1. Overview and results. We have assumed that the fork uj{a) admits a (con- 
tinuous and discrete) sum-of-poles representation H68|l . and the numerical investi- 
gations outhned in Sections 13.11 and 13.21 fcspeciallv the direct check in Section 
13.2.2(1 indicate that this assumption is valid. This sum-of-poles representation is 
quite similar in form to the one associated with integer-order Bessel functions and 
appropriate for wave propagation on a flat 2-1-1 dimensional spacetime. As such, 
we believe that Lemmas 3.4 and 3.5 in Ref. ^H] pertain to our representation 1(68(1 
of the Heun fdrk, demonstrating that for Hea > we may accurately approximate 
the kernel by a rational function, 

(99) ^" 

which is itself a sum of d simple poles. That these lemmas are indeed pertinent is 
an assumption. The approximating pole locations /3„ and strengths 7„ in 1(99(1 will 
be computed numerically. 

For our rational approximations -P(cr) and Q{<t) are polynomials with deg((5) = 
d = deg(P) -I- 1. Both P(cr) and Q{ij) depend on Z, pB, and k. As mentioned, we 
suppress this dependence throughout. Our choice P{a)/ Q{a) will be tailored to 
ensure that the relative suprcmum error 



(100) supj^gR 



|a)(a: + iy) - P{x + iy)/Q{x + iy) \ 
\Lj{x + iy)\ 



is smaller than a prescribed tolerance e. In practice, we set e to 10~^ or 10~^'^. 
Eq. 1(74(1 has motivated our interest in this error measure. We note that for the e 
considered, the second error term in (f75|l may be neglected so long as we aim for 
slightly better than the desired tolerance (for example, aiming for 7.5x10"" rather 
than 10~^°). As for the veracity of this statement, we offer the studies carried out 
in Sections I2.2.4I and 13 . 2 . 21 of the relative error |Aa>(i2/)|/|a)(iy)| in our knowledge 
of the true kernel. 

Throughout what follows, a; > is a fixed constant, and we have always set 
a; = in our work; however, for generality we retain x in the analysis. Our thinking 
here is that the kernel uj{ij) should be analytic in the righthalf cr-plane; hence Q{(t) 
should have zeros only in the lefthalf cr-plane (for simplicity, let us assume that 
P(cr) and Q{a) share no common zeros). If we set a; — 0, then the error above 
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Table 5. Number d of poles for given e. Number of poles 
needed to approximate the fdrk for = 15 with relative supre- 
mum error e. We find these numbers to be valid for both j = 0, 2. 
However, for the case j = 2 (gravitational radiation), one must 
assume I > 2. We have not yet examined the case .7=1 (electro- 
magnetic radiation). 



measures how well P{(j) / Q{a) approximates w(cr) along the Imcr axis, precisely the 
contour of integration appearing in the inverse Laplace transform. 

To find the desired rational approximation to the kernel, we solve a least-squares 
problem of the form 

2 



(101) minp^Q 



uj(x + iy) 



, Q{x + iy) 

that is we minimize the integral over the space of polynomials P(cr) and Q{(t) such 
that deg(Q) = d = deg(P) + 1. After we have numerically solved this problem and 
found P(cr) and (5(cr), we compute the relative supremum error ()1()0|I via numerical 
evaluations on a fine mesh. If this error is larger than the set tolerance e, then we 
increment d and try again. The upcoming Section 13.3.21 describes the algorithm 
we have used to solve (|101|l . We admit that for us the algorithm is half "black 
box." That is to say, although we shall describe it in some detail, we provide no 
hard analytical proof as to why the algorithm actually converges to a solution of 
the least-squares problem. Moreover, we do not discuss why minimization of this 
least-square error is associated with small relative supremum error. Nevertheless, 
we can check the accuracy of the final output and are able to achieve the desired 
tolerance. 

Before turning to a description of the algorithm, let us first give a quick overview 
of our results. In Table El we list the required number d of poles /?„ in H99() cor- 
responding to pb = 15 and two choices of the error tolerance e. We find these 
pole numbers to suffice for both the j = and j ~ 2 cases. Notice that d and I 
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Figure 32. Cut profiles for lowest-order kernels. Here 
we plot the (positive) Heun profile fo{x] 15) and the (negative) 
Bessel profile /(-i/2)(x; 15), both for pB = 15. To generate 
both we have used the one-component path method described in 
Section 12.1.31 The Bessel fork a;(_i/2)(a-; 15) associated with 
/(_i/2)(x; 15) is relevant for wave propagation in flat 2+1 dimen- 
sional spacetime. 



appear to grow at different rates, indicating increased performance in the large I 
limit. However, as we here consider only the relatively small bandwidth < Z < 64, 
we make this observation with some skepticism. We have extended the e = 10^^ 
side of the table and found that d = 14 is sufficient for I = 256. AGH have rigor- 
ously shown that for Bessel functions and flatspace wave propagation the number 
of approximating poles scales like 

(102) d~ 0(logzylog(l/£) +log2z/ + i/-Mog2(l/£)) 

as 1/ = Z + 1/2 ^ oo and e — > O"*". While it is beyond our mathematical ability to 
analytically establish this or similar growth for Heun functions, the growth indicated 
in Table IHl is not at odds with this result. Compression is also economical for low- 
order kernels which arc dominated by their continuous sectors. Tabic [S] lists d = 20 
as sufficient for e = 10~^° and 1 = 0. A numerical quadrature of the integral 
appearing in (|H5|l would require considerably more nodes to achieve this tolerance. 

As mentioned above, we have set x = in all cases (so that the least-squares 
problem is associated with the Imcr axis), and in particular we have done so for 
I = 0. This is perhaps remarkable in light of the fact that AGH found it necessary to 
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choose X nonzero while approximating the zero-order Bessel fdrk, a case for which 
the estimate (|74|l is not applicable. We stress that one need only deal with this 
zero-order Bcssel kernel (a circular rather than spherical boundary kernel) while 
considering wave propagation on flat 2+1 dimensional spacetime. To emphasize 
the difference between the I = Hcun fdrk and the i' = 1 + 1/2 = Bessel 
FDRK, we plot their cut profiles in FiG. |221 Note that for both cases the kernel, 
as represented by (|68|l . is determined solely by its cut profile. As documented by 
AGH in , loss of differentiability at the origin for the Bessel profile precludes a 
rational approximation of the zero-order Bessel fdrk by the methods described 
here, that is unless one chooses x positive and non-zero. For a non-zero positive 
X a more complicate estimate, valid only for finite time, must be used in place of 
(|74|l Fig. 1321 suggests that we may work with x = even for the lowest-order 
(Z = 0) Heun fdrk. 

3.3.2. Linear least-squares problem. To bypass the nonlinearity of the problem 
l)l()l|l . we switch to the linear iterative procedure 

2 



(103) 



lp(fc + l) Q(fc + 1) 



p{k+i)(x + iy) Q('^-+i)(a; + iy) , 

uj[x + ly) 



Q('=)(a; + i2/) Q('=)(2; + iy) 

where given Q'-'^^a) the task is to find the minimizing polynomials P^''~^^\a) and 
q('=+i)(o-), Xo commence the iteration, we need a degree d polynomial Q^°^(cr) 
as an initial guess, and we describe the construction of initial polynomials in the 
upcoming SECTION [3.3.31 

Assuming that we have Q^^^a) at our disposal, let us now discuss the solution 
to the least-squares problem p03() at the fcth level. To do so, we introduce the 
inner product 

/oo 
h{y)g{y)mk{y)dy , 
-oo 

where /i(y) and g{y) arc suitable functions and mk{y) ~ \ / \Q^'^'' [x + iy)\'^ is the fcth 
weight function. 

Claim: Set V{y) = P^^+^\x + \y), Q{y) = Q'^^+^'^ {x + \y) ., W{y) ^Lb{x + \y)., and 

( (x + i?/)"/2-i for n = 2, 4, 6,..., 2d 

(105) K{y) = I 

y (x + i2/)("^^'/2 W(2/) forri = l,3,5,...,2<i-l. 

Minimization of the integral l|103() is equivalent to 

(106) (-7'+QW,/i„)fc -0 
for n = 1, • • • , 2d. 

To verify the claim, first introduce the expansions V{y) = j=o "I" ^2/)"' 
2(2/) = Yl'j=o1j{^ + mV ^ assuming qa = 1. Next, consider a small perturbation 
of the mth coefficient Pm- Vanishment of the induced first-order variation of 
the integral 11U3|) may then be written as 



(107) {-V + QW,h 

2m+2)fc Spm + (— "P + QW, h2m+2) k^Pm — 0. 

Finally, exploit the freedom that Spm can be either real or imaginary, thereby 
reaching (|106(l for even n = 2m + 2. Eq. 1)1 06|) for odd n follows similarly upon 
introduction of Sqm perturbations. □ 
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With the claim (|106|l in mind, wc orthogonahze the 2d+ 1 functions 

(108) W{y), 1, {x + iy)W{y), (x + iy),..., 

{x + iyf-^W{y), [x + iyf-\ {x + iyYW{y) . 

The result is a family {gn{y) '■ n = 1, . . . ,2(i + 1} of orthogonal functions, with 
the last member g2d+i{v) the sought for solution ~'P{y) + Q(y)W(y) in (|106|l . We 
iteratively compute the functions in this family via the formulae 



(109) gn{y) 



{ W(y) 

1 - C2iW{y) 



{x + iy)gn-2{y) - 



inin{4,n — 1} 



c-njgn—j (y) 



for n = 1 

for 71 = 2 

for ?i = 3, 
. . . , 2d + 1 , 



where the real constants c„j are given by 



(110) 



for n = 2 and j = 1 



((gffn-2),gn-j)fc for n = 3, . . . , 2d + 1 
. {9n-3,gn-j)k and j = l,...,min{4,n-l} 



Here we use the notation {agn-2){y) for the function {x + \y)gn-2{y)- 

A few key observations confirm the correctness of these formulae. As suggested in 
FiGURE^of SECTiON i2.2.4l the profiles Ret2;(iy) and lu\Cj{iy) are of even and odd 
y-parity respectively, and for fixed a; > the profiles Recj(a; + \y) and Imd)(x + ly) 
also have the same y-parity. Therefore, along the path of integration used to define 
the inner product, W{y) has the form >V^(y) +iW~(y), where W"*" is even and W" 
is odd. Since x + iy is also of this form, it follows that all of the functions appearing 
in the list H1U8|I are as well. Now, the integration measure mk{y) appearing in 
the inner product {,)k is of even y-parity, showing that {,)k is real-valued on 
the space (I108|l . Whence the c„j are purely real. To see that c„j = for j > 4, 
first suppose that we have carried out the orthogonalization up to the function 
5„_2(y). In other words, we have g n-2{ y) orthogonal to gn^^iy), g^-iiy), . . . , gi(y). 
Owing to the pattern in the list H1U8|) . multiplication of gn-2{y) by x + iy yields 
a new function ((7(?„_2)(y) which is not in the span of gn~2{y) ign-ziy) , ■ ■ ■ ,gi{y)- 
Moreover, {agn-2){y) is automatically orthogonal to gn-j{y) for j > 4. Indeed, 
consider the identity 



(111) (x + iy)gn~2iy)gn-jiy) ^ 9n-2{y)[2xgn-j{y) -{x + iy)gn-j{y)] 

which implies that 

(112) ((CTg«-2),5«-j)fe = 2x{gn-2,gn-j)^. - {gn-2, (crg„-j))^ . 

Again owing to the pattern in the list (|108(l . the function {<ygn~j){y) = (x + 
^y)9n-j{y) is in the span of gn+2-j{y), gn+i-j{y), ■ ■ ■ ,5i(y)- Whence for j > 4 
we see that both terms on the RHS of the last equation vanish, giving Cnj = 0. 
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3.3.3. Numerical least-squares problem. We sketch how we have handled three is- 
sues which arise in numerically implementing the least-squares solve. 

Numerical quadrature. To handle the integration necessary to solve (|103|l via the 
described orthogonalization, we may change variables and write 

/oo /'7r/2 
-oo J-tt/2 

where h{y) represents any relevant integrand. This integral is numerically approx- 
imated as 

(114) ^ Ai/i(tan6l,)sec2 6l, , 

i=0 

where the 9i and A.^ are appropriate quadrature nodes and weig hts.^'^Due to the 
sum-of-poles form of the fork, all encountered integrands h{ta.n6) are periodic 
on [— 7r/2, 7r/2]. Moreover, all integrands are infinitely continuously differentiable, 
save possibly at 9 ~ Q. By analogy with integer-order Bessel functions, we expect 
regularity at = of order 21 + 1 ^H]. Whence for high I the trapezoid rule is 
effective and highly convergent. For low / the profiles Iieuj{iy) and lmLu{iy) are 
concentrated around and less regular at the y-origin. Therefore, for the first few I 
we employ an adaptive trapezoid quadrature rule, introducing more nodes around 
the origin. 

Evaluation of P{a), Q{<y), and their derivatives. Assume that we have produced 
a set of real constants {c„j : n = 2, . . . , 2d+l; 1 < j < 4} belonging to the fcth linear 
least-squares problem (|103f) described in Section 13.3.21 Further, assume that the 
solution polynomials pC^+^J (cr) and (5^^'+^) (a) are (close enough to) the polynomials 
P(cr) and Q{a) solving (|101|l . In other words, assume that the iteration in k has 
converged in some sense. In order to express the final approximating rational 
function as a pole sum H99|l . we must first know how to evaluate these polynomials 
and their derivatives at any a using only the Cnj constants.^'' This may be done 
via the following algorithm. Let Ai = 0, i?i = 1, A2 = —1, B2 = — C21, Ci = 0, 
Di = 0, C2 = 0, and D2 = 0. Next, for n = 3 until 2d + 1, set 

4 4 

(115) An — (jA„-2 — Cnj An- j , C„ = An-2 + O'Cn-2 — CnjCn-j , 

3 = 1 i=l 
4 4 

Bn = O-Bn-2 — ^ CnjBn-j , Dn = Bn-2 + ^^1-2 — ^ CnjDn-j ■ 

Then P{ct) = ^2^+1, Q{a) - Psd+i, P'{(t) = Cad+i, Q'(a) = Z?2d+i. With 
the ability to evaluate Q{<t) and Q'icr) in hand, we may use Newton's method to 
obtain the zeros /3„ of Q{a) which arc also the pole locations appearing in (|99|l . 
The corresponding pole strengths are then 7„ = P(/3„)/(5'(/3„). 



Do not confuse the quadrature weights A; with the fcth weight function mfe(j/) defining ( , }fe. 
In 11131 mj,{y) has been swept into h{y). 

17At each level in the iteration we must perform such an evaluation for Q^^\iy), thereby 
determining m}^{y) and in turn the fcth inner product. The relevant c„j come from the previous 
level. 
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Construction of the initial guess Q('')(cr). Recall that our sum-of-poles repre- 
sentation of the FORK features both a discrete sector and a continuous sector. With 
this in mind, we consider the following two cases. 

(i) d > Ni. The number of approximating poles is equal to or greater than 
the number of actual poles from the discrete sector. Therefore, one ap- 
proximating pole should correspond to each discrete-sector pole. Unless 
d = Ni, after each discrete-sector pole has been thusly represented, some 
approximating poles will be "left over." These should be placed on the 
negative Rccr axis in order to mimic poles from the continuous sector (cut 
contribution) . In case (i) we aim to naively mimic the actual pole locations 
found in the exact sum-of-poles representation. 

(ii) d < Ni. The number of approximating poles is strictly less than number 
actual poles in the discrete sector, and we do not attempt to naively mimic 
the pole locations found in the exact sum-of-poles representation. 

For case (i) we set Q'^^ct) = (ct — <;i){<7 — ';2) • ■ ■ (o' — ^d)- Then for all i < Ni 
we let Q = 15cri^i(15)/pB, that is we set q equal to the (numerically obtained) 
ith zero of Wi{15a] a) times the ratio 15/ps (we are also suppressing the I which 
could be on the cTi)- If there are any, the remaining d Ni locations cr^ are chosen 
as follows. For each relevant I (in practice for / < 10) we choose a fixed interval 
[a, h] C M<o lying on the negative Recr axis. This interval is "eyeballed" to ensure 
that [—6, —a] includes the essential support of the profile function /;(x; 15). We 
then send a 15a/ ps and b i-^ 15b/ pB- li d — Ni + 1, we let = (a + b)/2, and 
otherwise q — a + {i ~ Ni — l){b — a) / {d — Ni — 1) with i running from A^; -I- 1 to d. 

For case (ii) we set z = apB and then use the Besscl continued fraction expansion 
l)79|l . As mentioned there, Lenz's method may be used to evaluate this continued 
fraction. This method is described in the appendix of Rcf. ISI]. Via slight modi- 
fication it may also be used to return a partial denominator, which we use for the 
value (^(''^(iy). Indeed, even the most basic algorithm (see Ref. to compute 
continued fractions yields a partial denominator and may be used in this fashion. 

4. Evolution system and implementation of robc 

The goal of this section is to incorporate robc into the MacCormack algo- 
rithm, the finite-difference scheme we have used to simulate wave propagation on 
Schwarzschild blackholes. The first subsection describes the spacetime foliation and 
associated first-order system of evolution equations on which our numerical work is 
based. The second reviews the MacCormack algorithm in the context of this system 
of PDE, and it also addresses the issue of inner boundary conditions at the horizon 
H. The third subsection covers implementation of robc. The MacCormack algo- 
rithm is a simple, sccond-order-accuratc, predictor-corrector scheme, and wc show 
how ROBC fit nicely within the framework of the interior prediction and correction. 
In this final subsection we also touch upon some memory issues relevant to our 
implementation. 

4.1. Spacetime foliation and evolution system. 

4.1.1. Chosen foliation of spacetime into spacelike slices. Our numerical simulations 
of wave propagation on the Schwarzschild geometry feature physical coordinates. 
Moreover, for the time coordinate t and its associated foliation of the spacetime 
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into spacelike slices S we have chosen one closely related to ingoing Eddington- 
Finkelstein coordinates. The letter S is used both to denote the spacetime foliation 
and a generic slice of the foliation. There are several motivations for switching 
from the time T (introduced in SECTION II. 1.1(1 to t (defined in a moment). Besides 
removing the coordinate degeneracy present in the diagonal line-element (|14() . use 
of the time coordinate t also makes the inner boundary conditions at H particularly 
easy to set. Although we need not say more to motivate the upcoming discussion, 
further explanation may be found in Appendix C of 

Let us introduce the ingoing Eddington-Finkelstein coordinate system, here in 
terms of physical rather than dimensionless coordinates. We first pass to physi- 
cal coordinates by multiplying ((14(1 by an overall factor of 4m^ and subsequently 
sending 4m^ds^ i-> ds^. Next, we introduce the physical tortoise coordinate 

(116) = r + 2mlog [r/(2m) - l] , 
thereby reaching 

(117) ds^ ^ F{-dT^ + drl) + r^dd^ + sin^ dd<l)^ , 

where again F(r) = 1 — 2m/7'. Using the physical advanced time v — 2m.v = T + r*, 
we get the line-element in advanced Eddington-Finkelstein form, 

(118) ds^ ^ -Fdv^ + 2dudr + r'^dO'^ + sin^ Odcf? . 

Like the retarded time discussed in Section II.1.11 the advanced time v also labels 
characteristic surfaces, but ones which arc ingoing (cones which open up towards 
the past). From a numerical standpoint, we do not want to work directly with a 
null coordinate like v, rather a time coordinate T„ obeying w = T„ + r. Since by 
definition 

(119) V = (T + n -r) +r, 
the desired coordinate is 

(120) Ty^T + r^-r = T + 2m log [r/(2m) - l] , 

and wc refer to Ty ~ v — r as ingoing Eddington-Finkelstein time, although as 
indicated above it is the combination Ty + r which labels ingoing characteristics. 
Likewise, wc could introduce a time ^ u + r based on the physical retarded time 
u = 2m^ = T—r,,, but such a time variable would be less suitable for numerical work 
(see Appendix C of jj). Let lowercase t represent ingoing Eddington-Finkelstein 
time shifted by a constant as follows: 

(121) t^Ty- 2mlog [rB/(2m) - l] = r + 2mlog 

The family of level-t slices is the same as the family of level-T„ slices, and the 
combination t-\-r also labels ingoing characteristics. However, the t coordinate has 
been tailored to ensure that t = T at the outer boundary r = rs- 
By construction 



r — 2m 
7'B — 2m 



(122) 



dv ^ dTy + dr = dt + dr , 
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2 



2 • 2 

r sm 



sin^ 6'd^2 



and using this differential relationship, we rewrite H118() as follows: 

(123) ds^ = -Fdt^ + 2(1 - F)dMr + (2 - F)d7 

dt^ + (2 - F) L + ii— ^ 
The last line can be written 

(124) ds^ = -A^Mt^ + M2(dr + Vdtf + r^dO^ + sin^ 6d(t? , 
where the temporal lapse, radial lapse, and shift vector are the following: 

(125) N{r) = (1 + 2m/?')"^/^ , 

M{r) = (1 + 2m/?')^/^ , 
F''(r) = 2m/(2m + r) . 

In 127] York discusses the geometric interpretations of such variables in a general 
setting. For now, the main point is that the line-element (|124|l is nondegenerate at 
the horizon H , since N , M, and are all well-behaved at r = 2m. 

4.1.2. Wave equation with respect to chosen foliation. Let us derive the form of the 
wave equation in the Eddington-Finkelstein system {t, r, 9, 4>) of coordinates. This 
can be achieved via coordinate transformation at the level of pde, or by constructing 
the d'Alembertian associated with (|124|l . Choosing the first method, we start by 
rewriting H21|l in terms of physical coordinates, in order to reach 
-1 



(126) 1 



2m 

r 



QT2 



j,2 Qj, 



1 - 



2m 



dr 



2mj2^, 



0. 



In performing the near trivial coordinate transformation (r, p) = {t{T), p(r)), we 
have retained the symbol ipi to denote the new function satisfying the slightly 
new form (|126|l of the pde. Next, we introduce the function U{t,r) ~ ^'i{T{t),r), 
suppressing from here on out the subscript I which should be on U. Then, upon 
carrying out the coordinate transformation {T,r) = (T{t),r) on (|126|l . we find 

d^U 4m d^U ( 2m\ S^U 



(127) 



2m 

r 



dt'^ r dtdr 
2m dU 2{m-r)dU 
~~dt^ 



l{l 



2m 

r 

■f l)i7 



Qj.2 

2m fU 







Ot r^ dr r^ r-^ 

for the new form of the wave equation. Upon formal Laplace transformation on t, 
from this equation one recovers Leaver's normal form for the generalized spheroidal 
wave equation |44j . In order to trade this second-order equation for a first-order 
system of equations, we introduce an outgoing characteristic derivative of U . 

4.1.3. Outgoing characteristic derivative. Construction of the characteristic deriv- 
ative is based on an understanding of the relevant normal vector fields associated 
with the S foliation. Respectively, the following unit vector fields 



(128) 



e± = N-^ {d/dt - V'd/dr) , en = M'^d/dr 



are the future-pointing normal to the E foliation and the outward-pointing normal 
to concentric spheres within a given E slice j27j . These vectors are depicted in 
Fig. ESI 
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Figure 33. Concentric spheres in a E slice. The figure 
depicts two concentric splieres in a single Icvel-i slice as well as 
the vectors discussed in the text. The E normal e± points out of 
the slice into the future, while e\- points tangent to the slice but 
normal to the round sphere. These arrows may be extended as 
vector fields over the whole exterior region. 



With normal vector fields 1)128(1 . we introduce e+ = e±+e\- which is also depicted 
in Fig. 031 Then X = e+[U] is the characteristic derivative we seek, where [ ] 
denotes the standard operation of a vector field on a scalar function. More explicitly. 



1 dU 


f - 




dU 


N dt ^ 




N J 


dr 



(129) X 

We now take {U,X) as our basic variables. Since U and X belong to a single 
spherical-harmonic mode, both should carry I and m subscripts, but we continue 
to suppress these. 



4.1.4. First-order system of evolution equations. To facilitate numerical implemen- 
tation, we write (|127|l in a first -order form. The evolution equation for ?7, 

(130) dtU ^ iV' - M-^N)drU + NX , 

is just the definition of X. Eq. ((TTfjl now becomes the evolution equation for 

X. We write this equation in the general form 

(131) dtX =(y + M-^N)drX - N{X^' + 2X^)X + R^Xdra 

+ 2M-^NX^drU - R-^NI{1 + 1)U + 2mNR-^fU , 

where R is the areal radius, X^^ and X^ are certain characteristic variables of the 
background geometry, and a = N/{MR^) is the "dedensitized lapse." In this gen- 
eral notation one might denote our main characteristic variable 1)129(1 by X^ rather 
than plain X, but there is no need to do so here. In terms of the round-sphere area 
A, the variable R is defined geometrically as [A/(47r)]^/^, and, therefore, R = r for 
the Eddington-Finkelstein system. Moreover, in terms of the "Einstein-Christoffel 
variables" Km, Kr, Jm and fa written in Eqs. (2.2) and (2.3) of Ref. jSHl, we have 
X^^ = —Km — /m and X^ = —Kr + Jr. For the evolution at hand all of these 
variables are fixed, and in the numerical implementation they amount to function 
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calls given by 



(132) 
(133) 

(134) 

(135) 



i?(r) 
dra{r) 

X5(r) 



2fm + r) 



(2mr + 

r — 2m 
2mr + 



r2\2 



2m 



1/2 



2m^ — 5mr — 4r^ 
r(2m + r)2 



2m 



1/2 



We remark that the variable X'^ = e+ [log R] , whence describes the rate of change of 
the logarithm of areal radius along the outgoing null direction. Note that X:|^(2m) = 
0. Despite the somewhat general form (|131|) we have chosen for the wave equation, 
it is nothing but a tedious exercise in elementary algebra to substitute the given 
expressions for X, N, M, , R, X^ , X^, and dra into (|131|l in order to recover 
(iV times) Eq. p27|l . In performing this exercise, note right away that 



(136) 



yr 



M-^N = 1 



for the chosen coordinate system. Each of the evolution equations (|130I131|) has 
the form of a forced advection equation, and in particular we can write H131(l as 
follows: 



(137) 



dtX = drX + S{r,X,U,drU), 



where the source term S is defined by comparison of this equation with (|131|1 . Note 
that S also depends on the parameters I, j, and m. 



4.2. Interior of the computational domain and horizon: numerical details. 

4.2.1. MacCormack predictor- corrector scheme. Suppose that the radial mesh is 
the collection of Q + 1 nodes r^, where tq = 2m, rq = rs = 2iapB, and = 
2m + qAr. The spatial discretization step is Ar = [rs — 2m) /Q. Once the number 
Q of spatial subintervals is fixed, the temporal discretization step At is determined 
by fixing the Courant number At/Ar. Let the temporal mesh nodes be i„ = 
nAt, and define, for example, the mesh function [/^ = U{tn,rq). To simplify the 
presentation, our analysis holds J7^' as the exact value of U{t,r) at {tn,rq). Our 
task of numerical evolution on the interior of the computational domain is then to 
produce approximations, U^^^ and XJ^~^^, to the true values U^^^ and X^^^ (doing 
so for < q < Q), given in advance both and X^ (for < g < Q) as well as 
the exact geometry function calls Nq, Mq, Vq , Rq, [X^]q, [X^^]q, and [9ra]q listed 
in Eqs. ^7^^^^. Here, for example, Mq = M{rq) and [X^']q = X^'{rq). As 
these functions calls are time-independent, they need not carry the time superscript 
n. We choose the MacCormack scheme to accomplish our task. This scheme is in 
consistent conservation form as applied to a conservation law like the 1 + 1 scalar 
wave equation. Since our primary focus is numerical implementation of ROBC, it 
pays to keep the interior evolution simple in order to highlight the implementation. 



69 



Let us quickly review the MacCormack scheme in the context of the variables at 
hand. First consider the following predicted variables which stem from straightfor- 
ward differencing of H13UI137|I : 

(138) = U^' + At[{V^ - M-^N)^ D^Uj; + X;] , 

Here Dip denotes either a downwind or upwind difference stencil, 

(139) D^U- = {U^^ U-_,)/Ar , D+X^ = (x;^, X^) / Ar . 

To approximate the drU appearing in the source in the second equation of 
(|138|l . we use a centered difference stencil. The predicted variables Ug~^^ and 
are not second-order-accurate approximations to the values Ug '^'^ and X^^'^'^^ in 
question. However, the corrected variables 

(140) [7;+i = o.5{c7;+i + c/; + Atliv^ - m-^n)^ d+uj;+^ + x^+^j } , 

= o.5{x;+i + + At[D^xi'+^ + } , 

are indeed second-order-accurate approximations. To approximate the advcction 
derivative term in the correction phase (|140|l . wc have incorporated stencils D± 
opposite to those Dip used for the same term in the prediction phase (|138|l . The 
source term S^^^ at the next time step is built both with time-independent function 
calls and barred variables. Wc again use a centered stencil to handle the derivative 
term in the source, but now one 



(141) [drU]-+' = 0.5(c7;V+i' - U-+,')/Ar 
built with predicted values. 

4.2.2. Inner boundary conditions at H . From p31|l or H137|l we see that X propa- 
gates according to the cartoon H \ B, that is to say inward from B towards H, 
and at unit speed with respect to the coordinate time axis d/dt. Eqs. (|132I133I136|) 
then determine the following form for the X evolution equation at the horizon: 

(142) [d,X-drX]l^.^^^^ = 

- [NX^'X - r^Xdra + r-^Nl{l + 1)U - 2mNr-^j^U] \^,^^^^ , 

Notice the absence of the derivative term involving drU on the RHS. As X^{2m) = 
0, it has been dropped. Eq. H142|) shows that X propagates inward at H, whence 
there is no analytical issue of an inner boundary condition for X. 
The coordinate velocity appearing in H13U|) is 

(143) -V^ + M-^N ='—^, 
^ ^ r + 2m 

and it is monotonically increasing on [2m, oo) with < — y + AI^^N < l}^ This 
shows that U propagates according to the cartoon H B, at least for r > 2m. 



^^Based on 11431 . the unit-speed propagation of X , and the well— known stability properties 
of the MacCormack algorithm 1451 . we can therefore expect numerical evolution stability for 
At/Ar < 1. We have 1 on the RHS of this inequality, since 1, the speed of X everywhere, is larger 
than (pB — ^)/{pb the maximum coordinate speed for U over [2m, rs]. 
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More precisely, U propagates outward at speed |y — A/ ^A^| with respect to the 
coordinate time axis d/dt. Furthermore, at the horizon"'^^ 

(144) dtU\ , ^NX\ , , 

^ ^ I r— 2ni I r— 2m ' 

since — y + M~^N vanishes for r ~ 2m. Therefore, U propagates straight up the 
coordinate time axis at H, and there is no analytical issue of an inner boundary 
condition for U. Indeed, in tandem with the way 1)142(1 the variable X propagates, 
Eq. (|144|) implies that the region inside of r = 2m has no causal influence (insofar 
as the disturbances U and X are concerned) on the region outside of r = 2m. 
Heuristically, H acts as a one-way membrane. 

As for numerical inner boundary conditions at we extend the radial mesh past 
the horizon by adding a ghost node r_i = 2m — Ar. At both the prediction stage 
and correction stage, values at the ghost node are obtained via a simple copy-over 
from its righthand neighbor ro = 2m: C/'lj^ = C7j,'+\ X'lY = t^-l^ = 

and X'^^ = Xq^^. Of course, these are not the correct inner boundary conditions, 
but any resulting numerical error is trapped within the region r < 2m as shown 
by careful examination of the outlined MacCormack difference scheme at the point 
rp = 2m. Indeed, owing to the aforementioned identity [^]^]o — which kills the 
drll term in H131|) . the predicted variable ^q^^ depends neither on [/"j^ nor X^J^^. 
Likewise, since (y - M-^N)o = 0, the predicted variable also depends 

neither on U"i nor X"]^. Next, because X^^'l is a copy-over of X^+'^ (so that 
D^Xq^^ ~ holds), the corrected variable Xq^^ does not depend on mesh values 
at r_i. Therefore, we see that Uq~^^ also does not depend on values at r_i, in 
parallel with the fact that physically the region r < 2m cannot influence H. 

4.3. Implementation of robc. 

4.3.1. ROBC in physical coordinates. We may rewrite the convolution ((62|l in terms 
of the physical coordinates via division of the equation by an overall factor of 2m. 
The result is 

(145) X,(T,rB) + f![^'^^^ =rs'N{rB) ni{T ^ T';rB)MT' ,rs)dT' , 

M(rB)rB Jo 

where we have defined the physical tdrk 

u;i{T/{2my,rB/{2m)) 



(146) niiT;rB) 
Also appearing in H145|) is 

(147) X, = 



2ni 

19-0;, 1 di^i 



N dT M dr ' 

the characteristic derivative e+['(/';] of ipi along a null direction e+ which points 
outward and towards the future. 
We note that the null vector 

1 d 1 d 

^'''^ '--Ndf+Md-r 
points in the same direction as the one 

/i.nN 1 d f 1 V\ d 



^^In fact, the coordinate vector field d/dt, a Killing direction, is null at the horizon. 
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considered earlier. However, they have different extents and are related by a pure 
boost, 

(150) 8+ = e^e+ , 
through a hyperbolic angle 

(151) ^ = ilog ^^'^ '''^ 



1 - N-^MV 

Since the proper spatial velocity MV"^ /N = 2ni/r < 1 for r > 2m, the argument of 
the logarithm is positive on the exterior region. 

4.3.2. ROBC in Eddington-Finkelstein coordinates. The physical convolution H145|) 
is expressed with respect to spacelike slices which are level in the static time variable 
T. We now wish to express the convolution H145|) with respect to the spacelike slices 
S which are level in the ingoing Eddington-Finkelstcin time t. Since we have arrange 
for T = t at the outer boundary r — rs, we may express the physical ROBC given 
in (| 145(1 as follows: 

(152) e^"X(t,rB)+ ^ r^^N(rB) f n{t - t' ;rB)U{t' ,rB)dt' . 

M(rB)rB Jo 

To reach this equation, drop the I subscript on fi;, enact a trivial change of variables, 
and appeal to the equalities 

(153) Xi=e+[in]=e''e+[U] = e''X. 

The middle equality stems from p5()|l . and it of course holds at the outer boundary 
where we use the shorthand -Sb = ^^(t'b)- In this middle equality we have also 
dropped the subscript I. In (|152|1 we have chosen to retain N(rB) and M{rB), 
although they could be traded for N{rB) and M{rB) at the expense of introducing 
further boost factors. In any case, they are simply fixed constants, insofar as the 
ROBC are concerned. 

4.3.3. Approximate time-domain radiation kernel. From now until the end of this 
subsection we reinstate the policy of suppressing all factors of I. Moreover, for 
the most part we also suppress factors of ps, or as the case may be. The 
technique of kernel compression discussed in Section 13.31 yields an approximate 
FORK ^(cr) ~ P{a)/Q{a) of the form (j^H). The inverse Laplace transform £~^[^](t) 
of this approximate fork is an approximate tdrk 

d 

(154) e(r) = ^7«cxp(/3„r). 

Note that the d pole locations /3„ = (3^ + i/3^ and strengths 7„ = 7^ + 17^ are 
generally complex and the locations (are numerically observed to) lie in the lefthalf 
plane. In order to pass to a physical approximate tdrk S(t) ~ ^(t/(2m))/(2m), 
we set Pn = 2m6„ and 7„ = 2mc„, thereby reaching 

d 

(155) S(<) = ^ c„ cxp(6„t) . 

n=l 

We list a representative set of pole locations and strengths in Table IHl Note that 
these must be scaled for use in the physical kernel. For = 15 as in the table and 
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/3f + iPi 


-0.20807160 + iO 


+ 


-0.18971337+10.06900399 




-0.13304751 + 10.17772749 


/3f + 


-0.18971337 - 10.06900399 




-0.13304751 - 10.17772749 


7« + 


-0.07445512 + 10 


72" + 172 


-0.17152322 + 10.07630724 


I3 + 


-0.12569363 + 10.17804157 


i4 + hi 


-0.17152322 - 10.07630724 


l5 + 175 


-0.12569363- 10.17804157 



Table 6. Pole locations and strengths for a compressed 
KERNEL. Here the compressed kernel corresponds to j = 0, / = 4, 
PB = 15, and e = 10~^. Note that the number of poles is c? = 5 
In accord with Table El from Section |3.3.1I Both the in /3i and 
the In 71 correspond to output numbers from the compression 
algorithm which are of order 10^^^. 



a blackhole of mass m = 2, we would obtain the pole locations 6„ and strengths c„ 
corresponding to = 60 upon dividing the entries in the table by 4 = 2m. 

Of the d poles dsing will lie on the negative real axis and there will be dpair 
complex-conjugate pairs. So d = dsing + 2dpair always holds. For the compressed 
kernel shown in Table El we have dging = 1 and dpair = 2. We then break up 
as follows: 

(156) E{t) = Y,H,{t) + Y,G,{t). 

1=1 j=i 

In this decomposition 

(157) H,{t) = p, exp(K,i) , Gj{t) ^ 2 [mf cos(fcji) - rn} sln(fc}t)] exp(fcf i) , 

where Ki and /z^ respectively correspond to a pole location 6^+10 and strength c^+iO 
(each lying on the negative real axis), while A;^ + ikj and + Imj respectively 
correspond to a pole location + 1&^ and strength + Ic^ (each lying properly in 
the second quadrant). One should not confuse the fc^ + ifc^ with the zeros of the 
MacDonald function considered in Sections 12.1.31 and Ij.l.ll For the compressed 
kernel considered in TableElwe list these reordered locations and strengths in Table 
13 We also need to consider an auxiliary object, 

(158) Fj (t) = 2 [m^j cos(fcjt) + mf sln(fcji)] exp{kft) , 
not appearing directly in the kernel. We then use 

dsing „f rfpair 

(159) (S*t/)(t) = V / H,{t~t')U{t')dt' +y] / Gj{t^t')U{t')dt' , 

,^1 Jo ^-^1 Ja 

as an approximation to the convolution (fi * U){t) in (|152|1 . 
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2niKi 


-0.20807160 


2m(fcf + 


-0.18971337 + 10.06900399 


2ni(fc^' + iA:^) 


-0.13304751 + 10.17772749 


2ni/ii 


-0.07445512 


2ni(m{* + imj) 


-0.17152322 + iO.07630724 


2ni(m2^ + ) 


-0.12569363 + 10.17804157 



Table 7. Reordered pole locations and strengths for a 
COMPRESSED KERNEL. Hcrc wc show the reordered locations and 
strengths which render the compressed kernel In Tablel^lmanlfestly 
real. Compare these numbers with those In Table IHl A choice of 
mass must be made to determine these values. In this table and 
the text one should make a distinction between the mass (plain m) 
and the numerical pole strength rrij (Italic m and with a subscript). 



The update (.^ *U){t-\- At) plays a central role In our Implementation of robc. 
We write It as follows: 

(160) {E*U){t + At) ^ J E{t + At-t')U{t')dt' + E{t + At^t')Uit')dt', 

respectively referring to the first and second Integrals on the RHS as the history part 
h(S *U){t + At) and the local part l(5 *?/)(< + A<) of the updated convolution. In 
the same fashion, {Hi *U){t + At), {Fj *U){t + At), and {Gj *U){t + At) can each 
be split Into history and local parts. Computing the update (S * U){t + At) then 
amounts to computing the history and local parts of both [Gj * U){t + At) and 
{Hi * U){t + At). We discuss the local parts of these convolutions below, but here 
note that the history parts may be computed via the following exact Identities: 

(161) 

h(H, * U){t + At) = exp(K,Ai)(i?, * U)it) , 

h(Gj * U){t + At) = exp(/cf At)[cos(fc/At)(Gj * C/)(t) - sln(/s/At)(Fj * U){t)] , 

h(Fj * U){t + At) = exp(/cf At)[cos(fcjAt)(Fj * C/)(t) + sln(A:jAt)(Gj * U){t)] . 

These follow from elementary exponential and trigonometric Identities. The recur- 
sive evaluation of the history dependence afforded by identities H161|) leads to a 
considerable reduction in computational storage cost. A direct evaluation of the 
boundary convolution would require 0(t/At) memory locations to store the whole 
history of U along the timelike cylinder (see FiGUREnin the introduction). Our 
algorithm requires only the memory needed to store (the constituent pieces for) 
(S * J7)" as well as the numerical pole strengths and locations {2d = 2(ising + 4dpair 
real numbers). 

4.3.4. Incorporation o/ROBC into predictor-corrector algorithm. At t = t„ suppose 
that we are given (1) C/^" and X^" for < g < Q, (11) (S * J7)" = (S * J7)(t„; rs), and 
(iii) the necessary procedures to perform interior update as described in Section 
14.2.11 For the second assumption, we actually require Vi,_;/ that {Hi*U)^, {Gj*U)'^, 
and {Fj * [/)" are individually given. Again, for ease of presentation we assume that 
all of these present time-step expressions arc exact. Our implementation of ROBC 
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then amounts to the foUowing task: use (|152|l to produce second-order accurate 
approximations 

(162) U^+\ X'^+\ 

to the true values Uq+\ X'(^+\ and (S * Uy"+\ The notation (i7c/)"+i is short- 
hand for the pieces 

(163) {H~7u)''+\ (g7*^)"+S (ifrc/)"+i. 

We emphasize that Uq'^^ and Xq'^^ are the true values associated with a field freely 
radiating on a larger domain. 

We accomplish the task using the framework of the interior prediction-correction 
algorithm. Before entering the prediction phase, let us see what can be precomputed 
and stored for use during the course of the algorithm. First, we have 

fusing C^pair Casing C^pair 

(164) S(0)=^i/,(0)+^G,(0), E{At) = Y,H,{At) + J2G,{At) 

i—l J — 1 i— 1 i— 1 

at our disposal from t = Iq ~ onwards. Second, at the time step t = tn, we may 
appeal to the identities H161|) in order to exactly calculate and store the history 
parts uiHi * U)"+\ n{Gj * Uy+\ and n{Fj * C/)"+i Vi, j. With these we then 
exactly calculate the history part h(2 * C/)"+^ of the updated total convolution 

(s*;7)"+i. 

Let us now describe the prediction phase. Due to the way we have chosen the 
difference stencils in ()138|l . the boundary prediction Uq~^^ may be calculated from 
the same formula, 

(165) ^U^ + At[{V^ - M-'N)q {U^ U^-,)/Ar + X^] , 

used in the interior prediction. We then substitute Uq^^ into the approximation 
to (|152|l obtained by replacing with S, thereby obtaining the "predicted" char- 
acteristic variable at the next time step. Namely, 



(166) ^2+^ =r^^Nse-*« h(S * [7)"+^ + l(S * [/)"+i - C7^^+i 

where we use the shorthands Nb = N(rB) = l/M(rB) and 



(167) l(S * [/)"+i = 0.5At[E{At)W^ + E{0)U'^+^^ 



The last shorthand l(S * C/)"+^ is a prediction for the local part l(S * (7)"+^ of the 
updated total convolution (S * [/)"+^. This prediction is based on trapezoidal ap- 
proximation of the second integral appearing on the RHS of (|16()(l . In Section I5.1.2I 
we also consider what results from both rectangular and parabolic approximations 
as well as the trapezoidal approximation at hand. 
Turning to the correction phase, we first obtain 

(168) c7^+i=0.5{c75+i + C/5 

+ At[iV - M-'N)q {U^W - Ul+')/Ar + X'^+^J } . 

In this formula UqW is not available from the prediction phase, whence it is com- 
puted via the extrapolation 

(169) = 3C7^+i - Wlt\ + tJlt\ ■ 
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We discuss this detail of our implementation in the Appendix. Next, we compute 
the "corrected" variable 

(170) X^+i ^r^^Nse-''^ 
where as before 

(171) l(S^)"+i = 0.5At[E{At)Uj^ + S(0)C/^+i] . 
To end the correction phase, we similarly obtain Vi, j 

(172) l{F^)"+\ L{Gj7ur+\ L{H^r+\ 

and use these to update the constituent pieces of the convolution. For example, we 
take 

(173) (G^)"+i = h(G, * C/)"+i + l(g7^)"+' , 

as the desired approximation to the updated total convolution {Gj * U)"~^^. 

5. Numerical tests 

This section documents the results of several numerical tests of our implementa- 
tion of ROBC. Throughout, we consider a blackhole of mass m = 2 enclosed within 
an outer boundary of radius = 60 so that the horizon is located at 2m = 4 and 
Pb = 15. For this choice of the outer boundary B lies in the weak-field region 
(see the second paragraph of Section O. As mentioned in the introduction, we 
have chosen m = 2 only to have an example for which the mass is neither 1 nor 
i. Our numerical tests cover both j = 0, 2 cases as well as various values of /. 
Results carried out with pg = 20 are similar, but not reported. A more detailed 
investigation of these boundary conditions is forthcoming |67| . 

5.1. One— dimensional radial evolutions. Tests based on one-dimensional evo- 
lutions most effectively examine the issues of convergence and error for our ROBC, 
and this subsection is devoted to them. We first examine both issues in the con- 
text of a short-time numerical evolution. Subsequently, we further examine the 
issue of numerical error for our ROBC by considering two long-time evolutions. We 
also consider a third long-time evolution, although not in the context of an error 
analysis. 

5.1.1. Short-time evolution. Given fixed values for I and j, we have taken the initial 
data U{0, r) shown in the top plot of FiG. |M1 We have set [/(O, r) = /(r)/ r, where 



„(S *[/)"+! +l(S *[/)"+! -(7"+i 



(174) 



fir) = Aexp 



exp 



(r - a)2 {h- r)2 



for A = 50, a = 45, b = 55, and 5 = 20. The profile C/(0,r) is then compactly 
supported on [45, 55] and of C°° class. We choose X{Q, r) = — 1/(0, r)/r as the other 
initial condition. The resulting initial data set would be a pure outgoing were we 
considering / = and flat Minkowski spacetime. For the case at hand, whatever 
the fixed values I and j are, the initial data is not purely outgoing but will reach the 
outer boundary = 60 in a time (about 5) which is short relative to the crossing 
time of the domain [4, 60]. Recall that the spatial step size is Ar = (rs — 2m) /Q 
and the radial mesh points are rq = 2m -|- gAr. Now we have the particular values 
ro = 2m = 4 and rg = rg ~ 60. We arrange for a spatial discretization such that 
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Figure 34. Free evolution on a larger domain. The pa- 
rameters for the evolution are A^/Ar = 0.5, m = 2, ^ = 64, and 
J = 2. There are 32768 subintervals dividing [4, 116] and 16384 
subintervals dividing [4,60]. 



r2Q = 116. In particular, if ri6384 = 60, then wc have ^32768 ~ 116. Notice that 
112 = 116 — 4 is twice 56 = 60 — 4, and so [4, 116] is twice as large as [4, 60]. 

In order to generate an accurate reference solution, we numerically evolve the 
initial data on the larger domain [4, 116] until t = 38.5. For this evolution (and 
all others in this subsection) we choose a Courant factor At/Ar = 0.5, and at 
= 116 we use simple copy-over boundary conditions, since the disturbance 
does not reach this farther outer boundary in the short time t — 38.5. The result of 
our numerical evolution on [4, 116] is a numerical solution C7^™, described as "free" 
since it arises from a free evolution. Any numerical error in stems solely from 
the finite difference scheme and not from boundary conditions. The bottom plot 
in Fig. |^ depicts the result of such an evolution for the case / = 64 and j = 2 
(gravitational radiation). In both this plot and the plot of the initial data above 
it, we have indicated tb = 60 by a vertical line. The free solution that we compute 
corresponds to 2^^ subintervals of [4,60] and so 2^^ subintervals of [4,116]. Over 
the smaller domain [4,60], we have estimated that 

(175) ||A[/f™|U = sup{|[/f^<=<= - : q = 0, . . . , 16384} ~ 3.6 x 10"^ , 

where the f7*''"° are exact mesh-point values of the true solution. Notice that this 
supremum error measure is computed with truncated arrays (the full range of q is 
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Figure 35. Error of the solution [/Robc^ p^^^- ^j^g 
solute supremum error of the robc solution for various j, I, and Q. 
In all cases, these errors have been computed with the free solution 
corresponding to 2^^ subintcrvals of [4, 60] in place of the true so- 
lution. 



to 32768). We have obtained this estimate by a suitable convergence analysis of 
the truncated free solution and Richardson extrapolation |62| . 

Using our e = 10~^° ROBC, we may also evolve the initial data to t = 38.5 
directly on the smaller domain [4, 60], thereby generating another numerical solution 
jjROBC Furthermore, we may generate yet another numerical solution [/SOBc 
using Sommerfeld outer boundary conditions (SOBC) in place of our ROBC. By 
SOBC we mean the local boundary conditions defined by setting the numerical 
kernel Ei — throughout the evolution. The short-time numerical tests at hand 
then amount to numerical examinations of the convergence and accuracy of both 
jjROBC g^jjj [ySOBC relative to the free solution [/f''®'' (here serving in lieu of the 
true solution). Let us first consider the issue of convergence in the next paragraph, 
and afterwards turn to accuracy. 

Suppose fjROBC jg obtained for Q = 2^, with the supremum error |j A?7^°^'^||oo 
between [/ROBC g^j^j jjircc given by the formula 

(176) |lAC/"o^^|loo - sup{|C/™ - C/et^l : g = 0, . . . , 256} , 
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Figure 36. Comparison of robc and sobc. Here j = 2, I — 
64, and all profiles have been computed with 16384 subintervals of 
[4,60]. In the top plot the [/^O^C solution is the dashed line. 



where we note that 64 = 2^^/2^. In general, a stride of 2^^/2^ must be used to 
evaluate || AJ7^°^^||oo, when the array f7^°^^ corresponds to Q = 2'= for A: < 14 
(and U^'^'^'^ corresponds to 2^^ subintervals of [4, 60] as set up above). Let us switch 
to MATLAB notation for such norms, so that (|176|l becomes^" 

(177) II AJ7^°^^||oo = norm(i7^o^^(0 : 256) - C/f'™(0 : 64 : 16384), inf) . 
As we increase the number of subintervals used to compute {/ROBC^ .^g ^g^j^ check 
to see if the robc numerical solution initially converges to the free one. The 
plots in Fig. 1351 document this convergence for a small sampling oi ] and / values. 
Due to the equal scales of their horizontal and vertical axes, these plots indicate 
that our implementation of robc is second-order accurate as expected. Using L2 
errors instead of supremum errors, we obtain the same convergence rates as those 
indicated. Plots analogous to those in FiG. I^Hl indicate that for any choice of j 
and I values in our sampling the corresponding numerical solution [/SOBC jQ^g j-^ot 
converge to the free solution. 

Let us now consider the accuracy of both the robc and SOBC solutions, fixing 
for both a spatial resolution of 2^^ subintervals of [4, 60], the same resolution as for 
jjbee^ For the choices j = 2 and I = 64, FiG. 1311 depicts the numerical solutions 



We have slightly abused this notation, since MATLAB requires that the first entry of an array 
is labeled by 1. 
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Figure 37. Comparison of robc and sobc. Here j = 0, / = 0, 
and all profiles have been computed with 16384 subintervals of 
[4, 60]. In the top plot the [/SOBC solution is the dashed line. 



[/i^OBC aj^d C/SOBC ^3 ^^11 the errors AC/«obc ^ ^robc _ ^hoo ^jjSOBc ^ 
[/SOBC _ jjhec^ Fig. 123 depicts the same plots, but for j = and / = 0. For the 
sake of comparison, the top plots in these figures superimpose [/^^obc ^^^^^ jjSOBC^ 
That they differ is even evident to the eye, and a glance at the middle and bottom 
plots in the figures shows that the SOBC errors are many orders of magnitude larger 
than the ROBC errors. For the j — 2 and I = 64 case shown in FiG. 123 we note 
that II A[/^°^'^||oo (now computed with a stride of one) is about 2.38 x 10"''. This 
value is smaller than 3.6 x 10~^, our aforementioned estimate for the supremum 
error 

(178) ||At/f'-°'=||oo - norm([/'"'=(0 : 16384) - [/'™<=(0 : 16384), inf) 

in the free solution. In Table |H1 we list these error measures for the other values of 
J and / in our small sampling. In all cases the robc and free numerical solutions 
agree to within our confidence in the free solution. On the contrary, the SOBC and 
free numerical solutions do not agree to within this confidence, as is evident from 
the errors || A?7^'-'^"-'||oo also listed in the table. 

5.1.2. Long-time evolutions. We consider three long-time evolutions, taking for 
the first two j = 2, Z = 2, and At/Ar = 0.5 for all runs. Moreover, for the 
first two we choose initial data J7(0,r) = f{r)/r with /(r) as before in 1)174(1 . but 
now with for A ^ 10, a = 5, b ^ 15, and S = 20. The profile U{0,r) is then a 



I 


IIAC/'-'^lloo 






||A[/«^^^||oo 


64 


3.6 X 10-"^ 


2.46 X 10-' 


2.38 X 10- ' 


9.80 X 10-^ 


33 


9.7 X 10-' 


5.97 X 10-« 


5.65 X 10-*^ 


7.21 X 10-^ 





1.9 X 10-^ 


1.77 X lO-'J 


1.78 X lO-'' 


2.88 X 10-4 



Table 8. Error measures for short-time test. The table 
lists various supremum error measures for the j = case (the 
corresponding errors for j = 2 look the same). The first column 
of values for || AJ7^'''^''||oo are estimates for the error in the free 
solution, and we have obtained these estimates via a convergence 
analysis of the free evolution. In the last column, for example, the 
values II AC/^*^^^||oo measure the difference between the free and 
SOBC solutions when both are computed with a spatial resolution 
of 16384 subintervals. We hst two columns for || A;7^°^^||oo- The 
first corresponds to e = 10~^ and the second to e = 10~^° robc. 
We see essentially no difi^erence between these at the resolution of 
this test. 



compactly supported bump function of unit height on [5, 15], an interval overlapping 
the circular photon orbit at 3m = 6 |52|- We again set X(0,r) = —U{0,r)/r in 
order to crudely mimic outgoing initial data. Such data describes a wave packet 
that, although initially close to the horizon, somewhat escapes from the blackhole. 
For both evolutions, we consider a numerical test of our ROBC. Each one amounts 
to comparing a numerical solution JJ^OBC generated with our ROBC to a numerical 
solution generated on a larger domain, now cither four or sixteen times the 

size of [4,60]. 

These long-time tests show that our ROBC yield accurate numerical solutions 
corresponding to certain celebrated physical phenomena. In addition, the first 
long-time test further elucidates the error properties of the ROBC solution relative 
to the free one. Assuming that [/RO^C jjbee ^^^^ obtained at the same spatial 
resolution, any difi^erence ||t/^'-'^^ — C/^^^Hoo between them stems from three pos- 
sible sources: (i) use of an approximate kernel in the integral convolution at the 
boundary, (ii) the integration rule used to evaluate the integral convolution, and 
(iii) error in the boundary values of the solution U itself (which are fed into the 
convolution). Source (iii) has the potential for feedback. Our first long-time evolu- 
tion and experiment examines the competing influence of these sources. Based on 
this examination, we will conclude (to no surprise) that, given our second-order- 
accurate scheme for interior update, our trapezoidal implementation of the ROBC 
is optimal. 

First long-time evolution: quasinormal ringing. For the first test, we evolve the 
data on [4,228] until t ~ 201.25, thereby obtaining a new free solution. The top 
plot in Fig. 1381 depicts this free solution. Notice that the leading edge of the front 
has advanced from r = 15 to about r = 200, and that overall the disturbance falls 
off rapidly with decreasing r. However, the bottom plot, a graph of the base-10 
logarithm of the absolute value of rU^'^°° (with multiplication by r correcting for 
1 /r fall-off in amplitude) , exhibits the behavior called quasinormal ringing in the 
physics literature. This behavior is also depicted in FiG. 123 which shows the free 
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Figure 38. Free evolution on [4,228] at t = 201.25. 



evolution of the same data on [4, 900] until t = 700. As the leading front advances, 
it leaves exponentially decaying oscillations in its wake. The phenomenon is better 
described by examining the time history of the solution at a fixed radius, as is 
done in FiG. 1401 The complex frequency corresponding to this behavior is the 
least-damped (or fundamental) I = 2 quasinormal mode of the m = 2 blackhole 
in question. Kokkotas and Schmidt record a value of 0.37367 — iO. 08896 for the 
product of this frequency and the blackhole mass [HS]. Dividing this value by 2, 
we find 0.186835 — iO. 04448 as the least-damped / = 2 quasinormal mode of our 
black hole. The estimate tt/O. 186835 ~ 16.8 indeed corresponds to about half a 
wavelength in the Icfthand plot of FiG. 001 (as measured from one spike to the 
next), and the damping coefficient 0.04448/ In 10 ~ 0.019 matches minus the slope 
of the decaying amplitude in the same plot. 

To further examine the error properties of the robc solution, we generate the 
numerical solution [/ROBC [A,m] for t = 201.25, that wc may compare it with 
jjfrcc j-estricted to this smaller interval. The outer radius = 60 is marked by a 
black line in FiG.|2Hl Using matlab notation, for Q subdivisions of [4,60] we have 
the numerical solution [/^'-'^"-^(O : Q), and for AQ subdivisions of [4,228] we have 
the numerical solution J7^''°°(0 : 4Q). We may then compare the truncated solution 
U^'^^iO : Q) directly with U^'^^'^{Q : Q) by computing the supremum error between 
the two. Doing so for a sequence of spatial resolutions, we construct the column 
labeled trapezoidal in TableEl By changing the spatial resolution with /S.t/ Ar — 0.5 
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Q 


rectangular 


trapezoidal 


parabolic 


256 


1.6571 X IQ-s 


7.2426 X 10-^ 


6.8587 X 10-^ 


512 


8.3088 X 10-" 


1.7819 X 10- ^ 


1.6877 X 10- ' 


1024 


4.1663 X lO-'* 


4.3814 X 10-« 


4.1495 X 10-** 


2048 


2.0843 X 10-" 


1.0856 X 10-« 


1.0278 X 10-« 


4096 


1.0416 X 10-" 


2.7074 X 10--^ 


2.5627 X 10-** 


8192 


5.2053 X 10-^ 


6.7626 X 10-"^ 


6.4011 X 10-1" 


16384 


2.6023 X 10-^ 


1.6888 X 10-1" 


1.5985 X 10-1" 


32768 


1.3011 X 10-^ 


4.2149 X 10-11 


3.9891 X 10-11 



Table 9. Supremum errors for first long-time test. We 
list errors norm(U^'^^'^ {0 : Q) — ?7^''°°(0 : (5),inf) corresponding to 
(forward) rectangular, trapezoidal, and parabolic implementations 
of the ROBC. 



fixed, we also change the number of temporal steps taken in generating each solu- 
tion; however, we always compare ROBC and free solutions generated with the same 
number of temporal steps. Taking the base-2 logarithm of successive ratios of the 
trapezoidal errors in Table|5| we see a second-order convergence. Such convergence 
indicates that source (i) from the paragraph before last is not the dominant contri- 
bution to Ijc/ROBC _ ^frccii^^ Indeed, we have the same approximating numerical 
kernel regardless of mesh resolution. Therefore, for the numerical experiment at 
hand (which is based on our trapezoidal implementation of the ROBC) this error 
measure stems chiefly from the aforementioned sources (ii) and (iii). 

To further resolve the influence of these error sources, we have carried out the 
experiment using the rectangle rule (both forward and backward cases) , rather than 
the trapezoid rule as in (|167|l . to evaluate the local part of the boundary Laplace 
convolution. The corresponding errors for either case, with those for the forward 
case listed in the flrst column of Table El decrease at a first-order rate. Therefore, 
we conclude that for a rectangular implementation the dominant contribution to 
jjjjRQBC _ [/frooji^ source (ii). We have also carried out the experiment using 
Simpson's rule (parabolic rule) to approximate the local part of the convolution. 
Doing so involves a certain parabolic interpolation to get half-time-step numerical 
values for the field U. In this case, the corresponding errors, listed in the third 
column of Table again decrease at a second-order rate, although these errors 
are slightly smaller than their trapezoid counterparts. If (ii) were the dominant 
contribution to || [/RO^c _ ^ ^^j, ^^j, trapezoidal implementation of the ROBC, 

then switching to a parabolic implementation should lead to error decrease at a 
fourth-order rate. Since it does not, we conclude that our trapezoidal ROBC are 
not the dominant source of error (assuming a second-order interior scheme) . These 
considerations suggest that there is little point in improving the integration rule 
used for the Laplace convolution, unless one also implements a higher order scheme 
for the interior. 

Second long-time evolution: decay tails. Our second long-time, one-dimensional, 
test evolution shows that the ROBC yield a reasonably accurate numerical solution 
even for very late times relative to the crossing time of the domain. After the 
quasinormal ringing dies down, the solution at a fixed radius r is known to behave 
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Figure 39. Free evolution on [4,900] at t = 700. 

as order ^-(2^+3)^ provided the initial field configuration is not time-symmetric 
|68j . and this behavior is seen in the righthand plot of FicEHl Heuristically, this 
late-time decay tail stems from the back-scattering of outgoing waves off of the 
background curvature. FicOHl depicts the solution at < = 700, and the bottom plot 
shows that by this late time a sizable region has entered the power-law decay regime. 
Note that rs = 60 is again marked by a black line in the figure. For Q subdivisions 
of [4,60], there are 16Q subdivisions of [4,900]. Therefore, we truncate the full 
numerical solution C/^™(0 : 16Q), and like before compare [/^™(0 : Q) directly 
with [/^'^^'^(O : Q). Doing so for a sequence of spatial resolutions, we construct 
Table ^1 Due to round-off error, these errors do not decrease at a second-order 
rate. However, with the spatial resolution 4096 indicated in the table's fifth row as 
an example, we note that the supremum relative error^^ 

(179) norm((C/f''='=(0 : 4096) - [/^^^^(O : 4096))./C/f^'=<=(0 : 4096), inf) 

is the reasonably accurate value 2.1019 x 10""*. For SOBC, rather than 6.0051 x 10~^^ 
for the absolute error and 2.1019 x 10~^ for the relative error, we get the numbers 
2.9569 X 10-9 and 5.2668 x 10^. Fig. BTl depicts the solution at t = 700 on [4,60]. 
The top plot superimposes both C/f''='=(0 : 4096) and U^'^^'^{0 : 4096) to show that 
they match to the eye. 

In the bottom plot of FiG.|iT]we show a modified U^°^'^{0 : 4096). This profile 
has been generated by "turning off" the cut contribution to the kernel. More 



The ./ is MATLAB notation for component-by-componcnt division. 
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Figure 40. Quasinormal ringing and decay tail. For the 
described second long-time evolution, we record the time history of 
the field at the fixed position = 60 (the vertical line in FiG. I39|) . 
The lefthand plot depicts quasinormal ringing. The righthand plot, 
a zoom~in of the late-time history of the field, features logj^Q t to 
highlight the onset of the power-law decay. Each plot is a 
history of the numerical boundary value J7q^^^. 



precisely, to obtain this profile, in 1159|) we have by hand set Hi{t — t') = 0, Vi 
throughout the evolution. For this I = 2 case, the kernel S{t) = Gi{t) in (|156|l then 
corresponds to a single conjugate pair of poles. From this experiment, we might 
infer that the phenomenon of decay tails is handled by the cut contribution to the 
kernel. As seen in Sections 13.11 and 13.21 (but there for j = 0) , this contribution 
stems from a continuous distribution of poles. 

Third long-time evolution: decay tails revisited. Here we consider an experi- 
ment described by Allen, Buckmiller, Burko, and Price [HEIj in order to further 
demonstrate that our robc capture the phenomenon of decay tails. In terms of 
the standard retarded time m = T — we introduce a function f{u) = \u{u — 
8m)/(16m2)]8 for < u < 8m, with f{u) = otherwise. Since we work with 
Eddington-Finkclstcin coordinates, whereas |3f)j worked with static time coordi- 
nates, we set u{t, r) = t + r — 2r* — C . Here C is a constant, chosen to ensure that 
C + 2mlog[-l + C/(2m)] = 0, and in turn u(t, C) = t. [Before in (fHT)! and the 



norm 


^^Kuau^O : 256) - (0 : 256), inf) 


1.8166 X 10"^'^ 


norm 


[/HUh!U(o : 512) - C/*rcc(o ; 512), inf) 


1.1427 X IQ-^^ 


norm 


[/HUBU(o : 1024) (0 


1024), inf) 


7.6272 X lO-i*^ 


norm 


[/HUBU(o : 2048) (0 


2048), inf) 


6.2889 X 10-1** 


norm 


[/HUBU(o : 4096) (0 


4096), inf) 


6.0051 X IQ-^*^ 


norm 


[/HUBU(o : 8192) - [/'"<=(0 


8192), inf) 


5.9483 X 10-1** 



Table 10. Supremum errors for second long-time test. 
Here we are using matlab notation for the supremum error as 
explained in the text. 



analysis foUowing that equation, t denoted a time variable which differs from this 
one by a constant, but no matter.] We obtain an initial j,l ~ scalar wave packet 
as U{0,r) = f{u{0,r))/r, and also complete the data accordingly with X{0,r) ob- 
tained from H129(l . That is, we get X{0,r) by computing e+[f{u{t,r))/r] and then 
setting t = 0. Therefore, the data describes an essentially outgoing pulse which at 
t = is of unit height and roughly supported on [2m, C]. Next, on a numerical 
domain [2m, 30m], we evolve the data until t = 1000m with our e = 10~^° ROBC, 
along the way recording the time history of the product rC7^°^*-^ of radius and nu- 
merical field at the fixed radius C (more precisely, at the radial mesh point closest 
to C). Here we multiply by radius (whereas in FiG. I40l we did not) in order to have 
better agreement with [SS]. The resulting history is shown in FiG. 021 We stress 
that our domain size is more than 16 times smaller than the one considered in |36| . 

We believe this experiment shows that our ROBC indeed resolve decay tails. This 
is not to say that the results of [SEI are incorrect. Those results are applicable for 
the type of boundary conditions considered, namely approximate conditions which 
are not history dependent. We intend to carry out other more careful tests of our 
ROBC along the lines spelled out in [SHI . Results will be presented in a forthcoming 
work |67j . 

5.2. Three— dimensional evolution. Our final test is three-dimensional, and 
considers the scenario of a wave-packet striking the computational boundary at 
a shallow angle (or at least not on a perpendicular). Such a scenario is known to 
be difficult test case for radiation boundary conditions in flatspace. We carry out 
this test only for the scalar .7 = case. 

5.2.1. Description and set~up. Consider the bump function 
(180) 

/(t,a;,y,z) = Aexp ( ^ ) / exp 




|r — vt — ro| 6)2 (6 — |r — vt — Fq |) 
where 

(181) |r - vt - ro| = [{x - - xof + {y - Vyt - yof + {z - v^t - zof]^^^ . 

Were we considering wave propagation on flat Minkowski spacetime, / would be 
an exact solution. Were this the case, the bump function would have center Tq = 
{xo,yo, zq) at t = and would travel without distortion and at unit speed in the 
direction v = (vx ,Vy,Vz). Here we are assuming that v'^ + Vy + f ^ = 1 . 
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Figure 41. U^'"" and u^obc ^ ^ 700. The modified f/ROBC 
solution has been obtained by removing the cut contribution to the 
TDRK as described in the text. 

Let us interpret (x, y, z) as the "Cartesian coordinates" belonging to the Edding- 
ton-Finkelstein coordinates (r, 9, (f) from Section 14.1.11 (the same spatial coordi- 
nates as in Section II .1.11) . That is to say, x = rsin^cosc/), and so on.^^ Doing 
so gives us a difFeomorphism between the initial Schwarzschild time slice S and 
three-dimensional Euclidean space E'^, itself viewed as the initial surface of flat 
spacetime. On E"^ we obtain the data /|t=o and i9t/|t=o which is then taken over 
as initial data on S. The Schwarzschild evolution of this data is not the spacetime 
function H18U|) . Nevertheless, on physical grounds, one expects this evolution to be 
qualitatively similar to the flat~spacetime evolution (|180|l of the same data, pro- 
vided that (i) the center (iq, 2/0, zq) is of large enough radius [xq -I- 2/0 + ^oY^"^ ^^'^ 
(ii) the support b relative to that center is small enough. 

We have written a three-dimensional spectral code to evolve such data. The ba- 
sic idea, implemented numerically, is to obtain a data set, t/;m(0,r) and X;m(0,r), 
for each spherical mode {I, m), a set we then evolve via the one-dimensional purely 
radial scheme described in Section 0] (although throughout that section we sup- 
pressed the Im subscripts on these mode variables). The exact data for each mode 



'These are not the well-known Schwarzschild system of isotropic coordinates 1521 . 
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Figure 42. Long-time history of the field rU^'^^^ at r = 
C ~ 2.5569m. Here we show the resuh of the third long-time 
experiment with j, ^ = and initial data inspired by |3f)j . C is the 
radial value corresponding to vanishing tortoise coordinate. More 
precisely, here we consider the history at the radial mesh point 
~ 2.5572m closest to C. This plot may be compared with Figure 
2 in 121]. The dash-dot line is the base-10 logarithm of 4.57(m/t)'^, 
similar to the matching curve considered in . This plot exhibits 
the power-law decay of the field. We stress that this plot has been 
obtained with the outer boundary radius tb = 30m. In carrying 
out the experiment, we set m = 2. 



are obtained by harmonic analysis as follows: 

(182) t/,,„(0,r)= / / f{Q,x{r,9,^),y{r,9,^),z{r,e,cj)))Yi.rn{O,^)sin'0d9dcj), 

Jo Jo 

and 

(183) Xi^iO, r) = / / g{0, x{r, 0, cj)),y{r, 9, cj)),z{r, 9, cb))Yi„,{9, 0) sin^ 9d9d(f>, 

Jo Jo 

where g = dtf+r^^{xdxf+ydyf+zdzf)- After all one-dimensional evolutions have 
been carried out, we have a collection of Uim{t, r) modes (here viewed as exact for 



rsine -100 ^ .^qo 



Figure 43. Equatorial cross section of C/^™'' at t = 38.5. 
The initial data for this evolution is described in the text. The 
outer, middle, and inner circles are r = 116, 60, and 4. For this 
plot (Q, I, J) = (2048, 33, 32), so that the shown spherical domain 
(determined by 4 < r < 116) corresponds to a (4096,33,32) dis- 
cretization. The oscillations spread out over the whole angular 
domain result from the harmonic truncation of the initial data. 

ease of presentation) with which we can construct the three dimensional solution 

oo / 

(184) U{t, r,e,c^) = Yl ^'™(*' ^)Yhn{0, 0) 

/=0 m=-l 

via harmonic synthesis. Numerically, we obtain f\t=o and g\t=o on a spherical-polar 
grid 

(185) {(r,, 0„ 0,) : < g < Q, 1 < i < /, 1 < j < J} 

by exact function calls. For each we then perform a harmonic analysis on the 
restriction of the initial data to the two-dimensional spherical grid {9i,(j)j)- We 
have used Adams and Swarztrauber's SPHEREPACK 3.0 to numerically perform 
all harmonic analysis and synthesis, and in particular those subroutines (shagci, 
shagc, shsgci, and shsgc) from SPHEREPACK which take the 9i as Gaussian points 
in colatitude. 

We work with an angular grid corresponding to (/, J) = (33, 32), with the first 
number odd so that the equator is a Gaussian point. The spectral truncation in 
SPHEREPACK is triangular, and with this angular resolution we only sample har- 
monic modes for I < 32. Prior to the harmonic analysis described above, we 
perform a preliminary analysis followed immediately by synthesis. This initial pro- 
cedure corresponds to truncation of the full harmonic series for f\t=o and g\t=Oi a-nd 
it does modify of the described initial data. The resulting data is now spread out 
over the whole angular domain, although it remains concentrated at (a:o,yoj^o)- 

5.2.2. Short-time three-dimensional evolution. The three-dimensional numerical 
test we now describe is somewhat analogous to the one-dimensional short-time test 



Q 




||AC/«'^^^|l2 






256 


1.7593 X 10-'' 


4.9594 X 10-*^ 


7.3416 X 10-^ 


8.3218 X 10"^ 


512 


1.6646 X lO"'' 


1.1387 X 10"" 


7.3137 X 10-^ 


8.2895 X 10"^ 


1024 


3.4185 X 10-« 


2.7872 X 10-^ 


7.3131 X 10-^ 


8.2724 X 10-^ 


2048 


7.9344 X 10- ^ 


6.8970 X 10-« 


7.3131 X 10-^ 


8.2656 X 10-^ 



Table 11. Error measures for the three-dimensional test. 



described in Section 15. 1.11 We consider a spherical domain bounded internally by 
r = 4 and externally by r = 116. On this large domain, we choose initial data as 
described in the last subsubsection. That is to say, in addition to the discretization 
sizes (Q, / = 33, J = 32) of our three-dimensional domain, in H180|) we also choose 
the vectors (xq, j/Oj -zo) S'lid {vx^Vy^Vz) as well as the values of A, 6, and b. We 
then numerically evolve this data until t = 38.5. In setting up the initial data, 
due care must be taken to ensure that (i) the support of the data is contained 
within r < 60 and (ii) in the short-time t = 38.5 the propagating solution spills 
over r = 60 but does not reach r = 116. To satisfy these criteria, we have chosen 
a bump function (|180|l centered at {xo,yo, zq) = (—40, 20, 14), with initial velocity 
{vx,Vy,Vz) = (—0.87,-0.4,-0.3). Moreover, we have chosen A = 1, 6 = 12, and 
5 = 0.1. 

The byproduct of our evolution is a three-dimensional free solution 1/^'^'^° at 
t — 38.5. As before, any numerical error in 1/^'^°'^ stems solely from the finite- 
difference scheme and not the boundary conditions. Having generated such a free 
solution from the initial data described above, we show the result in FiG. 031 a 
depiction of the equatorial cross section of U^'^^'^. As mentioned, the number / = 33 
of di points ensures that 9n = 7r/2, allowing for such a figure. Our initial data has 
been designed to yield an evolution with an interesting cross section. 

The numerical test consists of again evolving the data, but now on the smaller 
spherical domain determined by 4 < r < 60 and using the e = 10^^° robc. 
As a result, we get another three-dimensional array [/f^o^c -^^rjjjgjj -^^g compare 
with a suitably truncated U^'^'^'^. Note that now both [/RO^c g^^^^ jjirco ^^.^ three- 
dimensional arrays of physical values, whereas in the first subsection they stood 
for one-dimensional arrays of spectral values. Wc consider two error measures 
corresponding to the difference /^jj^obc _ [|R0BC _ [/free^ 'pj^g gj.g^ jg gij^piy ^j^g 

absolute supremum error 

(186) II At/i^OBC||^ ^ ^MK?r - ^iTi I : < g < Q, 1 < z < /, 1 < J < J} , 

while the second is a discretization of the L2 error 
(187) 

.7 I Q 



]AC/"OBC|, . - 1 



vol(I]) 



1/2 

E E E - u!-::,\'M,ri sln^.ArA^.A,/, 



where the radial lapse factor Alq = (1 + A/rgY^"^ ensures that this discretization 
corresponds to proper integration. The volume of the computational domain is 
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(188) vol(I]) = 4:71 (1 + A/ry/^r^dr ~ 9.48255 x 10^ 
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Q 


ROBC 


SOBC 


256 


50.00 


46.46 


512 


49.45 


48.15 


1024 


43.21 


41.85 


2048 


19.55 


18.18 



Table 12. Relative times for three-dimensional evolu- 
tions. Here we list the relative run times as percentages for the 
various Q values. With Q = 1024, for example, we generated the 
free solution on the larger domain via a run of 1817.26s. The corre- 
sponding ROBC and SOBC runs on the smaller domain took 785.18s 
and 760.61s respectively. The entry 43.21 in the table corresponds 
to 100 X 785.18/1817.26, and so on. 



Since the 0-mesh is comprised of Gaussian points 9i, we must also include the 
appropriate Gaussian weights A^i in the sum (both the 9i and A0i are obtained 
with the routine gaqd from SPHEREPACk). We emphasize that unlike the situation 
in the first subsection where all error measures have been computed in the spectral 
space of spherical harmonics, these error measures are computed in physical space. 

In Table ITTI we list errors for the scenario shown in FiG.031 doing so for radial 
discretizations corresponding to Q = 256, 512, 1024 and 2048 (so that the free 
solution has been computed with radial discretizations corresponding to 2Q = 512, 
1024, 2048 and 4096). In the table we have also listed the corresponding Sommerfeld 
errors, which have been computed as in ifT^ and but with U^'^'^'^ in place 

of [/ROBC^ Notice that the ROBC give rise to perfect second-order convergence, 
whereas the SOBC do not. Moreover, the ROBC numerical solutions are several 
orders of magnitude more accurate than their SOBC counterparts. It is also of 
interest to consider the relative computational cost of both ROBC and SOBC. As 
SOBC are purely local boundary conditions, they are of course less expensive than 
our ROBC. However, despite being strikingly more accurate, we find our ROBC to 
be nearly as cheap as SOBC, as is evident by the relative timings listed in Table [T^ 

6. DISCUSSION 

We start this final section by comparing our results with those of our main 
reference |18j . thereby gaining some perspective on what we have achieved. After 
that comparison we discuss possible applications and extensions. 

6.1. Comparison with the work of Alpert, Greengard, and Hagstrom. Our 

rapid implementation of ROBC on the Schwarzschild geometry has been inspired by 
Alpert, Greengard, and Hagstrom's implementation on flatspace. Although we have 
followed their approach closely, our implementation falls somewhat short of theirs. 
Let us point out why this is the case by comparing both implementations from 
theoretical and numerical standpoints. 

From a theoretical perspective our work is less satisfactory than that of AGH, 
chiefly because — as mentioned in the introduction — we are unable to give a rig- 
orous asymptotic analysis of our implementation. For each angular / index AGH 
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deal with a Bessel fork which is theoretically known to admit a sum-of-poles rep- 
resentation, one corresponding to a purely discrete sum in the case of 3+1 wave 
propagation and to a discrete and continuous sum in the 2+1 case. Moreover, due 
to the exhaustively studied properties of Bessel functions, they also start with a 
wealth of useful information about how the pole locations, pole strengths, and (in 
the 2+1 case) cut profile for the representation behave in all conceivable asymptotic 
regimes (which include certain scaling properties as / becomes large). With tight 
control over where the physical poles accumulate and the behavior of cut profile, 
they are able to borrow ideas from the fast-multipole method in order to replace 
(discrete and continuous) physical pole sums with (fully discrete) approximate pole 
sums of fewer terms. As a result, they rigorously prove that the sum-of-poles rep- 
resentation for the Bessel fork admits a rational approximation, one uniformly 
valid in the righthalf frequency plane and exhibiting exponential convergence as 
the number of approximating poles is increased. For the scenario we have consid- 
ered, such an analysis would seem out of the question. Indeed, although wc have 
provided extremely convincing numerical evidence that the Heun fdrk admits a 
sum-of-poles representation, one strikingly similar to the 2+1 Bessel case, even this 
is conjecture from a theoretical standpoint. The battery of asymptotics needed to 
theoretically prove that such a representation admits a rational approximation in 
the style of AGH is certainly well beyond the author's knowledge of special functions. 

Our numerical implementation of robc also falls somewhat short of AGh's, in 
that they considered a higher bandwidth (their I < 1024 compared to our I < 64) 
and a smaller best error tolerance (their e = 10^^^ compared to our e = 10^^"). 
However, we believe that on this count our work is closer to being on the same 
footing with theirs. Indeed, insofar as numerical implementation of flatspacc ROBC 
is concerned, their elegant asymptotic analysis proving the existence of rational 
approximations is somewhat beside the point. Ultimately, their compression algo- 
rithm (which yields the desired rational approximations) relies only on the ability 
to evaluate the Bessel fdrk along the imaginary axis of the frequency plane. For 
Bessel functions, which obey certain order recursion relations, such evaluation can 
be efficiently done via the continued fraction expression H79() following from such 
relations. We stress that such evaluation requires no knowledge of the sum-of- 
poles representation for the kernel. Although we have no such continued fraction 
expression with which to evaluate the Heun fdrk, we have seen that our integra- 
tion method is almost as accurate (in the sense spelled out by Section I2.2.4|) . For 
Bessel kernels we have observed that the numerical path integration required by our 
method is more expensive than continued fraction evaluation, and all the more so 
as the order ^ + 1/2 gets large, although wc have not made a systematic comparison 
of the two methods. However, the cost of evaluation is almost beside the point, 
since in principle any extra cost associated with our method need only be incurred 
once. Of true importance is accuracy, and through the use of extended precision we 
believe it possible to build Heun kernels corresponding to e = 10~^^ and through a 
bandwidth of 1024.^'^ We point out that, insofar as both gravitational wave astron- 
omy [HHl and the post-Newtonian approximation |64l 17(1] of the gravitational field 
are concerned, I values well below 64 are the ones primarily relevant to gravitational 
wave observation. 



Of course, had AGH taken advantage of extended precision with the continued fraction 
method, they presumably could have pushed beyond even their reported numbers. 
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Another numerical shortcoming of this work is that we have mainly considered 
certain discrete values of the outer boundary radius ps, whereas AGH allow for 
any value. We hope to extend our results to a continuous interval in pB, say 
[15,25], via construction of our rational approximations at Chcbyshev nodes in 
\/pb and subsequent interpolation. The only reason we have not done this so 
far is that the straight compression algorithm described in Section 13.31 yields 
rational approximations which do not appear to vary that smoothly with l/ps 
(it is a nonlinear minimization after all). We anticipate that this problem can be 
mitigated. Indeed, using the described compression, one could construct a rational 
approximation at a Chcbyshev node £}g = l/p]^ near 1/15, and as a result obtain a 
number of pole locations [the Q in HlQl|l ]. For the other nodes, one could simply fix 
a Q by suitable scaling of the collected locations for and then solve a straight 
least squares problem for P in (|101(l . This straight least -square solve would be the 
same problem as (llUlll but with Q fixed. In any case, we do plan to build up some 
form of our ROBC which is valid over a continuous ps interval. 

6.2. Potential applications and extensions. We briefly describe several poten- 
tial applications of our ROBC. The first arises in the theory of non-spherical stellar 
collapse j52| . A massive star, having expended its nuclear fuel, will collapse into 
a blackhole, with gravitational waves radiated in the process. Provided that the 
initial configuration of the collapsing system is not highly non-spherical, the evolu- 
tion of such waves are governed by the Regge- Wheeler and Zerilli equations. Our 
ROBC are therefore applicable to such scenarios. The classic work on perturbed 
Oppenheimer-Snyder collapse is that of Cunningham, Moncricf, and Price |71ll72j . 
but we hope to carry out further numerical investigations with our ROBC. Another 
field of recent interest where our ROBC should find application is the theory of stellar 
perturbations, in particular oscillation modes of stars as sources for gravitational 
radiation (see for example |73l 1741 17^ V As a third application, we anticipate that 
our ROBC can by applied to the Cauchy-perturbative matching scheme of Rupright, 
Abrahams, and RczzoUa a numerical approach to ROBC for full general rel- 

ativity. The idea is to match the Cauchy evolution of the full Einstein equations 
to a set of one-dimensional, linear, radial evolutions of perturbative modes on 
Schwarzschild, with the matching taking place at an extraction two-surface large 
enough to ensure the validity of perturbation theory. In this approach, each of 
the modes is separately evolved on a large radial domain until the next time-step, 
allowing for both wave-form extraction and updating of the interior solution at 
the extraction surface. The equation governing the evolution of such modes is 
related to the one we have considered. Rupright et al employ Sommerfcld outer 
boundary conditions for their radial evolutions. We expect that our ROBC could 
improve their method, perhaps even doing away with the separate radial evolutions 
altogether (the idea here would be to impose our ROBC directly at the extraction 
surface, provided it indeed sits well within the perturbative region). 

Let us mention several other arenas where the issue of ROBC might be examined 
along the lines laid down here. The most obvious would be ROBC for time-domain 



They are easily able to allow for any value for the following reason: The flatspace radial wave 
equation arising from Laplace and spherical— harmonic transformation (essentially the modified 
Bessel equation) can be expressed solely in terms of sr, the product of Laplace frequency and 
radius. 
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wave propagation on more complicated blackholes: the Reissner-Nordstr0m so- 
lution (charged Schwarzschild), the rotating Kerr solution, and the charged Kerr 
solution (also known as Kerr-Ncwman) While tackling the problem of robc 
for a metric as complicated as the Kerr solution would be very difficult, one might 
first obtain ROBC for an approximate solution (to the Einstein equations) which 
exhibits rotation, such as the Brill-Cohen metric (see p. 699 of [^21 and references 
therein). On a different track, robc might be investigated for other types of pde 
on the Schwarzschild geometry, such as curved-space versions of the Klein-Gordon 
or Euler equations. Implementation of robc for fluid equations on blackholes 
would have application in realistic studies of blackhole accretion. The methods 
of this work could be extended to the Schrodinger equation with a potential, like 
the Coulomb potential, which is not of compact support. (In his dissertation 
Jiang studied ROBC for the Schrodinger equation via the AGH approach, but only 
considered compactly supported potentials.) Implementation of robc in such a 
setting would involve the special function theory of Coulomb wave functions or 
the related Whittaker functions (both incarnations of confluent hypergeometric 
functions) |58l 1591 IHOj . Finally, we remark on a more practical possibility, namely, 
extension of these methods to the ordinary wave equation but expressed in prolate 
or oblate spheroidal coordinates 00] . Ref. [70j addresses computational issues asso- 
ciated with evaluating spheroidal wave functions. The issue of robc in this setting 
would be not unlike the issue for Kerr blackholes. In any case, the spheroidal wave 
equation is essentially the confluent Heun equation, so one might expect our meth- 
ods to be applicable. Radiative boundary conditions which are tied to spheroidal 
surfaces might prove useful in modeling phenomena with slender geometries, such 
as those arising in antenna design (see Ref. j77) and references therein). 

Yet another extension of our work would involve coupling robc to high-order 
interior schemes, such as the one described in |78| and |79j . but without reduction of 
accuracy. This is a delicate issue pertinent to the simulation of waves on flatspace as 
well as on blackholes. We have presented a second-order accurate implementation 
of ROBC, but it would be beneficial to achieve a fourth or higher-order implementa- 
tion. Even with the exact robc in hand, it is a nontrivial problem to couple them 
to a high-order interior scheme, as a naive coupling gives rise to numerical bound- 
ary layers which — upon spatial differencing — spoil the order of accuracy. Indeed, 
such problems occur when trying to couple exact Drichlet or Neumann boundary 
conditions to well-known high-order schemes, such as Runga-Kutta schemes, and 
they have a long history j45)-|51j. For achieving a successful high-order implemen- 
tation of robc, a strategy based upon Picard integral deferred correction would 
seem promising (see |8()| as well as comments made by Minion |81p. 
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Appendix A. Modified MacCormack scheme 

In Section 14.3. 41 we use the extrapolation 1)169(1 in lieu of the predicted variable 
UqXi which is not yet available, and this is a salient feature of our implementation 
of ROBC. We stress that this extrapolation takes place only at the outermost mesh 
point tq = rB- Nevertheless, in order to investigate its validity, we ask whether 
or not it may be extended to a valid numerical scheme over the whole interior 
computational domain. 

For simplicity and in order to make precise statements about the scheme to be 
considered, let us focus attention on the simplest conservation law 

with a > 0, as a model of Eqs. ((13()|1 and 1)143(1 . Now consider the following modified 
MacCormack scheme. The prediction phase is unchanged. Namely, we take 

(190) = ^i{U- U-_,) , 

with /i = aAt/ Ar (not to be confused with the retarded time coordinate). However, 
for the corrected variable we set 

(191) C/^+i = [C7^+1 + - /.(2t/^+i - 3U^+l + U^+i)] . 

The three terms within the round parenthesis correspond to (C/^^i — C^g^^) in the 
straight MacCormack scheme. However, we use the extrapolation 

(192) 3C7;'+^ - Slf^+l + U'^q+l 

in place of U^^l- This modified scheme is second-order accurate, as can be verified 
by a calculation based on Taylor series. Moreover, it can be written as 

(193) ^ [/; - ^,[F{u-, c/,"_i, c/^2) - f{u-_,,u-_,, u-_,)\ , 

in terms of the flux function F{U, V, W) determined by 

(194) 2F{U,V,W) =3U -V - ti{2U ~3V + W) . 

Since F(U, U, U) = J7, the scheme is in consistent conservation form as applied to a 
conservation law like the 1-1-1 flatspace advection equation 03] . Following Minion's 
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description |H1] of a standard stability analysis, we have written a short matlab 
code (given below) which demonstrates that the modified scheme is stable for the 
time-step constraint aAt/Ar < 2/3 (or a number shockingly close to 2/3). The 
output plots for this code are depicted in FiGS. 1441 and 1^ The second plot shows 
that instability in the scheme stems from the highest modes. 

For the straight interior MacCormack scheme as applied to the system H13()|l and 
(|13H) of evolution equations, we have remarked in Section 14.2.21 that we expect 
evolution stability for At/ Ar < 1. Recall that this is the expectation because the 
variable X propagates everywhere with speed 1, while {ps — ^)/{pb + 1) < 1 is the 
maximum speed of U over the domain [2m, r^], as is evident from (|143|1 . We now 
assume that our stability analysis for the modified MacCormack scheme — carried 
out using the model equation (|189|l — also pertains to (|13()(l . Since the modified 
scheme is only used at the outermost boundary point, we let a = (ps — l)/(ps + !)• 
Then for pB ~ 15, we get a = 14/16, and in turn the constraint At/ Ar < 16/21 ~ 
0.7619 associated with the modified scheme. Since the rhs of this constraint is less 
than 1, the issue at hand is whether it is 16/21 or 1 which limits evolution stability. 
From numerical experiments, such as those in Section [3 we find that stability 
is ensured so long as At/ Ar < 1. Although the modified MacCormack time-step 
constraint is more restrictive, it is not a limitation for our numerical evolutions, 
presumably because the modified scheme is implemented only at a single point. 

"/. 

S R Lau 

Applied Mathematics Group 
y, Department of Mathematics 
University of North Carolina 
Chapel Hill, NC 27599-3250 
USA 

"/. 

26 October 2003 

"/. 

y, Matlab code: Modif iedMacCormack.m 
"/. 

Code produces two plots, each of which elucidates the time-step 
stability associated with applying a modified MacCormack scheme 

y, to the 1+1 linear advection equation, 

7. 

I U_t + a U_r = 0, 
"/. 

y, where a > is a constant velocity. Our method of analysis 
y follows M. L. Minion, "On the Stability of Godunov-Projection 
y Methods for Incompressible Flow," J. Comp. Phys . {\bf 123}, 435 
y (1996) . 

y 

y The modified MacCormack scheme is as follows. In terms of the 
y Courant-Friedrich-Levy factor 

y 

y mu = a Dt/Dr, 

y 

y we obtain the predicted variable 

y 

y barU(n+l,q) = U(n,q) - mu [U(n,q) - U(n,q-1)] . (A) 
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7. 

y, Next , we construct the corrected variable 
"/. 

tilU(n+l,q) = 0.5-[barU(n+l,q) + U(n,q) - mu 
•/. [2barU(n+l,q) - 3barU(n+l ,q-l) + barU(n+l ,q-2)] > . (B) 

7. 

7, The three terms within the square parenthesis correspond 

7o to [barU(n+l ,q+l) - barU(n+l,q)] in the straight MacCormack 

7o scheme. However, here we use the extrapolation 
7. 

7 3barU(n+l,q) - 3barU(n+l ,q-l) + barU(n+l,q-2) 
7 

7 in place of barU(n+l ,q+l) . This scheme is inspired by how 

7 ROBC have been implemented within the MacCormack scheme in 

7 S R Lau, "Rapid Evaluation of Radiation Boundary Kernels for 

7 Time-Domain Wave Propagation on Blackholes . " UNC-Chapel Hill 

7 Ph.D. dissertation (2003). See the end of Chapter 4. 

7 

7 To perform a standard Von Neumann stability analysis, suppose 

7 that U is periodic on [0,2pi], and represented by the discrete 

7 Fourier series 
7 

7 Q/2-1 

7 U(j) = Sum hatU(p) exp(j p Dr) , (C) 

7 p = -Q/2 

7 

7 where Dr = 2pi/Q and hatU is the Fourier transform of U. Then 

7 the range of xi = p Dr lies in [-pi, pi]. To define Amp(xi,mu), 

7 the amplification factor, we substitute the spectral form (C) 

7 into equations (A) and (B) above, in order to get the symbol 

7 Amp for one step of the method: 
7 

7 hatU(n+l,p) = Amp(xi,mu) hatU(n,p) . 
7 

7 Having given the relevant background, let us turn to the code. 
7 

7 step 1. Build up a two-by-two mesh grid with CFL factor mu 

7 along X-axis and angle xi = p Dr along y-axis . The argument 

7 of each trigonometric function which makes up the amplification 

7 factor Amp involves p Dr = 2 pi p/Q, so to test Amp we sample 

7 xi from -pi to pi at regular intervals (continuity of Amp 

7 ensures that we need not sample every xi) . 

7 

7 Choose discretization sizes for mu and xi arrays, and build the 

7 mesh grid. 

7 

Dmu = 0.001; 
Dxi = pi/16; 

[xi,mu] = meshgrid(-pi :Dxi :pi , :Dmu: 1) ; 
7 

7 step 2. Build corresponding two-by-two array of values for 

7 the Amp symbol particular to the modified MacCormack scheme. 
7 
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Expl = exp(-i*xi) ; 

Exp2 = exp(-2*i*xi) ; 

Exp3 = exp(-3*i*xi) ; 

terml = 3-4*Expl+Exp2 ; 

term2 = 2-5*Expl+4*Exp2-Exp3; 

Amp = 1-0 . 5*mu. *terml+0 . 5*mu. *mu. *term2 ; 

"/. 

y. step 3. Make the first plot. For each point in the mu array, 
y, we compute the maximum of |Amp(xi,mu)| over all xi . Note mu = 
y, always gives an amplification factor of unity, so maxAmpMod > 1 
y, marks instability. 

y. 

[J K] = size (Amp) ; 
for j = 1:J 

maxAmpMod(j) = max(abs (Amp( j , : ) ) ) ; 
end 

figure (1) 
hold off 

plot (mud : J, 1) ,maxAmpMod) 
hold on 

xlabel('\mu = a Dt/Dr') 

ylabeK 'max\-[Amp(\xi , \mu) : -\pi \leq \xi \leq \pi\}-') 
title ('Time — step constraint for modified MacCormack scheme') 
plot ((2/3)*ones(size( [0.8:0.1:3] )) , [0 . 8 : . 1 : 3] , 'k: ' ) 
axis tight 

y. 

y, step 4. Make the second plot. Here we plot contour lines of 
y, the modulus of the Amplification factor in order to see which 
y xi are the most amplified. 
1 

figure (2) 
hold off 

[CS H] = contour (mu,xi ,abs (Amp) ) ; 
xlabel('\mu = a Dt/Dr') 
ylabeK '\xi') 

title (' Contour lines of I Amp(\xi , \mu) I ' ) 

clabel(CS,H, 'manual') 

1 

y Press "return" to exit manual labeling. 

y 
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Time-Step constraint for modified MacCormack scheme 



3 



2.5- 



E 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

]X = a Dt/Dr 

Figure 44. First output from the given matlab code. 
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Contour lines of IAmp(^,)a,)l 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

^ = a Dt/Dr 

Figure 45. Second output from the given matlab code. 
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Appendix B. Numerical Tables for Compressed Kernels 

We list the pole locations /3„ and strengths 7„ for a few compressed kernels used in the 
numerical tests of Section |^ A table entry of indicates an output number from the 
compression algorithm which is less than 10"^*^ in absolute value. For those /3„ and 7„ 
appearing in conjugate pairs, we list only the one with positive imaginary part. 



n 




Pi 


n 


7^ 


7^ 


1 


-0.407381913421500E+00 





1 


-0.560438492794723E-06 





2 


-0.292310786462780E+00 





2 


-0.149307319758660E-04 





3 


-0.212438334101443E+00 





3 


-0.100427245970089E-03 





4 


-0.154351731119108E+00 





4 


-0.294347741100713E-03 





5 


-0.111621522695175E+00 





5 


-0.482077412607607E-03 





6 


-0.801606769182339E-01 





6 


-0.512454401482300E-03 





7 


-0.570590618998134E-01 





7 


-0.391967225466562E-03 





8 


-0.401964994812843E-01 





8 


-0.233213890995425E-03 





9 


-0.280020468146328E-01 





9 


-0.114815532845651E-03 





10 


-0.192813689223556E-01 





10 


-0.490292206150649E-04 





11 


-0.131170606343391E-01 





11 


-0.187703451634142E-04 





12 


-0.881053984705117E-02 





12 


-0.658176657111368E-05 





13 


-0.583815623973656E-02 





13 


-0.214158551560558E-05 





14 


-0.381311155357965E-02 





14 


-0.651562115228849E-06 





15 


-0.245325710136527E-02 





15 


-0.186157771876132E-06 





16 


-0. 155499824057403E-02 





16 


-0.501212404555003E-07 





17 


-0.970762744563191E-03 





17 


-0.128840007496407E-07 





18 


-0.585846647089220E-03 





18 


-0.314418019702290E-08 





19 


-0.373076210020314E--03 





19 


-0.569152368353698E-09 





20 


-0.183147313598350E-03 





20 


-0.244872495848780E-09 






Table 13. Kernel for j = 0, ^ = 0, ps = 15, e = 10" 
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n 




PL 


1 


-0.347849789139467E+01 





2 


-0.328013420888283E+01 


0.122459742396181E+01 


3 


-0.277834690894067E+01 


0.221807916777226E+01 


4 


-0.216780239221729E+01 


0.288987484597755E+01 


5 


-0.160611143055369E+01 


0.328338093803925E+01 


6 


-0. 1 16571976808040E+01 


0.348766466680969E+01 


7 


-0.870203833982946E+00 


0.359360642350047E+01 


8 


-0.654574374170026E+00 


0.373250677729775E+01 


9 


-0.383775970479126E+00 


0.391836743018780E+01 


n 


In 




1 


-0.351882332982489E+02 





2 


-0.288247787077508E+02 


0.145963809805392E+02 


3 


-0. 159408675098397E+02 


0.197009045186916E+02 


4 


-0.579079784769916E+01 


0.165503367808947E+02 


5 


-0.888910445944541E+00 


0.110514097610705E+02 


6 


+0.617612362898097E+00 


0.649132105745073E+01 


7 


+0.488858645232752E-01 


0.363187756402738E+01 


8 


-0.604551383223953E+00 


0.358904868249248E+01 


9 


-0.352475685403885E+00 


0.379107013210578E+01 



Table 14. Kernel for j = 2,1 = 64, pB = 15, e = 10"^". Complex 
conjugation of entries 2-9 gives entries 10-17. 



n 




Pi 


1 


-0.375616176922446E+00 





2 


-0.252285897920276E+00 





3 


-0.171458781191188E+00 





4 


-0.116562490243471E+00 





5 


-0.764331990276028E-01 





6 


-0.468091057981411E-01 





7 


--0.263730137927373E-01 





8 


-0.125652994567027E-01 





9 


-0.947795178946719E-01 


0.599312024947409E-01 


n 


In 




1 


-0.942815440763951E-05 





2 


-0.366046310052021E-03 





3 


-0.374027383588918E-02 





4 


-0.872734265927278E-02 





5 


-0.147189136342128E-02 





6 


-0.501356988668342E-04 





7 


-0.973423621068051E-06 





8 


-0.728807025058112E-08 





9 


-0.894836172990982E-01 


0.620643548936884E-01 



Table 15. Kernel for j = 2 J = 2, pB = 15, e = lO^^". Complex 
conjugation of entry 9 gives entry 10. 
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